MM for Molecular Dynamics Simulations: 6

Subscriber access provided by Iowa State University | Library

Dynamics

Adaptive Partitioning QM/MM for Molecular Dynamics Simulations: 6. Proton Transport through a Biological Channel Adam W. Duster, Christina M. Garza, Baris O. Aydintug, Mikias B. Negussie, and Hai Lin J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b01128 • Publication Date (Web): 14 Jan 2019 Downloaded from http://pubs.acs.org on January 21, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Dec 27, 2018 Revised for J. Chem. Theory Comput.

Adaptive Partitioning QM/MM for Molecular Dynamics Simulations: 6. Proton Transport through a Biological Channel Adam W. Duster, Christina M. Garza, Baris O. Aydintug, Mikias B. Negussie, and Hai Lin* Chemistry Department, CB 194, University of Colorado, Denver, Colorado, 80217, United States

Received Date: *Email: [email protected]

KEYWORDS: Adaptive partitioning, QM/MM, CLC, proton transport, water wire

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 43

ABSTRACT

Adaptive quantum-mechanics/molecular-mechanics (QM/MM) dynamics simulations feature on-thefly reclassification of atoms as QM or MM continuously and smoothly as trajectories are propagated. This allows one to use small, mobile QM subsystems, the contents of which are dynamically updated as needed. In this work, we report the first adaptive QM/MM simulations of H+ transfer through a biological channel, in particular, the protein EcCLC, a chloride channel (CLC) Cl–/H+ antiporter derived from E. coli. To this end, the H+ indicator previously formulated for approximating the location of an excess H+ in bulk water was extended to include Cl– ions and carboxyl groups as H+ donors/acceptors. Furthermore, when setting up buffer groups, a new “sushi-roll” scheme was employed to group multiple water molecules, ions, and titratable residues along the one-dimensional channel for adaptive partitions. Our simulations reveal that the H+ relay path, which consists of water molecules in the pore, a bound Cl– ion at the central binding site (Cl–cen) of the protein, and the external gating residue E148, exhibits certain mobility within the channel. A two-stage journey of H+ migration was observed: The H+ moves towards Cl–cen and is then shared between Cl–cen and nearby water molecules in the first stage and departs from Cl–cen via nearly-concerted transfer to protonate E148 in the second stage. Most of the simulated trajectories show the bound Cl– ion in the channel to be transiently protonated, a possibility that was previously suggested by experiments and computations. Comparisons with conventional QM/MM simulations revealed that both adaptive and conventional treatments yield similar qualitative pictures. This work demonstrates the feasibility of adaptive QM/MM in the simulations of H+ migration through biological channels.


2

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction The past decade has witnessed a rapid growth in the development of adaptive algorithms1-27 in the framework of combined quantum-mechanics/molecular-mechanics (QM/MM).10,

27-52

Adaptive

QM/MM reclassifies atoms as QM or MM on-the-fly continuously and smoothly during dynamics simulations, which is not allowed in conventional QM/MM. Adaptive QM/MM is particularly suitable for modeling ion solvation and transport, for which the ion and its immediate surroundings are often of primary interest and should ideally be described at the QM level of theory. With adaptive QM/MM, one employs a small, mobile QM subsystem (also called the active-zone) that follows the ion wherever the ion goes, and updates its contents dynamically as needed. The advantages are apparent: one does not have to decide which atoms will be QM or MM prior to starting the simulation, and a small QM subsystem facilitates the utilization of higher-level QM theory and/or longer simulation times, which may potentially lead to new insights. Recently, we extended the adaptive algorithms to treat a solvated H+ in bulk water.15,

18

Unlike

ordinary ions (e.g. Na+ and Cl–), a hydrated H+ is usually viewed as a delocalized charge or geometric defect. According to the Grotthuss mechanism,53-56 H+ transport is characterized as the propagation of the charge or geometric defect by reorganizing the network of covalent and hydrogen bonds rather than as the migration of the given H atom. This leads to a constantly changing identity of the transferred H+. For example, in bulk water, an excess H+ is generally represented by the hydronium O in an Eigen structure and by the shared H in a Zundel structure. The ever-changing identity of the moving H+ is a challenge to adaptive QM/MM, where the center of the QM subsystem is typically located at the ion of interest. To solve this problem, we have proposed an H+ indicator as a means for keeping track of the hopping H+.15 The H+ indicator’s location is a function of the coordinates of the donor and its surrounding acceptors, and varies smoothly as the H+’s solvation structure changes between the Eigen and Zundel structures. We have demonstrated that the geometry-based H+ indicator is a faithful


3


Page 4 of 43

representation of the hydrated H+, yielding dynamics and potential of mean force in test calculations15, 21, 57

that are comparable to those58-64 obtained with another commonly used H+ representation, the

(modified) center of excess charge.59, 62, 63, 65 The H+ indicator is essential to the success of the adaptive QM/MM simulations of a hydrated H+ in water, because the mobile QM zone can be centered on the indicator.15, 21, 57 In the present contribution, we report the first adaptive QM/MM simulation of H+ transport through a biological channel. Compared with the simulations of an excess H+ in water, biological channels have two distinct features: (i) H+ relay through a channel is basically a movement along one-dimensional hydrogen-bonding network as opposed to the three-dimensional movement in bulk water, and (ii) titratable ions and functional groups may be part of the H+ translocation path (known as the H+ wire). Accordingly, two new treatments are introduced in this work. First, we describe a “sushi-roll” strategy that increases the efficiency of adaptive calculations by dividing the channel into groups of crosssectional slices. Second, we generalize the definition of the H+ indicator for an excess H+ in bulk water to include monoatomic ions and the titratable carboxyl group. These two new developments are important for future adaptive QM/MM studies of H+ transport through biological channels. The biological channel studied here is the E. coli CLC transport protein (EcCLC).66-72 The chloride channel’s (CLC) genes are widely expressed in from humans to prokaryotes,66,

73

and CLC transport

proteins form a ubiquitous group of Cl– ion channels and Cl–/H+ antiporters. In human body, CLC proteins are associated with many critical physiological and cellular processes such as aiding extreme acid response, cell-volume regulation, and muscle contraction.66, 74-76 As an antiporter, EcCLC regulates influx of Cl– and efflux of H+ in a 2.2:1 ratio across membranes (Fig. 1).68-70, 77-79 Unlike conventional antiporters, where substrates take turns occupying a single pathway, H+ and Cl– can simultaneously bind to EcCLC, and their transport pathways overlap.80,

81

Intriguingly, mutation of a single residue, the

external gating E148, can convert this antiporter to a Cl– channel.67,

68

The operating mechanism of


4

Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


EcCLC is very complicated, which is to a large extent due to the coupling of the Cl– transport and H+ relay in opposite directions. Many experimental67, 68, 70, 72, 76, 79-96 and computational9, 61, 80, 93, 94, 97-128 efforts have been devoted towards elucidating the acting mechanisms of this transporter. For example, it has been suggested by experiments70 and by computations115 that the Cl– ion bound in the central binding site (Scen) may be transiently protonated by the migrating H+. Notably, a recent theoretical study127 based on the concept of Markov State Model129-132 employed both computational and experimental data to optimize rate coefficients and revealed that the Cl–:H+ exchange ratio of 2.2:1 can arise from the combination of the kinetics for the possible transitions between different states.

Figure 1. Illustration of the pathways for Cl– (green) and H+ (red) transport, bifurcated at E148. The membrane is represented by blue continuum. The two subunits of EcCLC (PDB code 1OTS67), which can function independently, are displayed as helixes in yellow (chain A) and gray (chain B), respectively. In chain A, the Cl– ions are displayed as green spheres in the central (closer to E148) and intracellular binding sites (closer to the intracellular solution). Critical residues are highlighted in balls and sticks: S107 and Y445 for Cl– transport, E203 (internal gate) for H+ transport, and E148 (external gate) for both. The inlet gives a detailed view of the pathways with water molecules in the pore (positions taken from MD simulations93), where only the Cl– ion at the central binding site (Cl–cen) is shown. Color code: C, cyan; H, white; N, blue; and O, red.

In this work, we study one specific scenario of H+ transfer: with one Cl– ion present in the central binding site (Scen) and with the internal gating residue E203 protonated, one excess H+ hops from one


5


Page 6 of 43

water molecule near the internal pore entrance to the deprotonated external gating residue E148 through water molecules in the pore.109 In this particular scenario, the EcCLC protein functions like a passive channel and serves as a “toy” model for the demonstration of our adaptive QM/MM algorithms. While the relatively short path for H+ relay (~15 Å across) is not ideal for demonstrating the time-saving advantages of adaptive QM/MM over conventional QM/MM, it does allow us to perform reference conventional QM/MM simulations that would otherwise be too expensive to validate the adaptive QM/MM calculations. We focus on two aspects of the H+ relay dynamics: the dynamic pathways that the H+ travels through and the interactions between the H+ and Cl–cen. This paper is organized as follows: The methodology will be described in Section 2, followed by the computational details in Section 3. The results will be presented and discussed in Section 3, and conclusions will be drawn in Section 4.


6

Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2. Methodology 2.1. Permuted Adaptive Partitioning Scheme with Electrostatic Embedding In this work, we employ the modified permuted adaptive-partitioning (mPAP) scheme,9 which is a variant of the permuted adaptive-partitioning (PAP) scheme.3, 6 The family of the PAP algorithms has demonstrated superior numerical stability in MD simulations.3, 6, 9, 15, 27 Details about these schemes can be found in the original literature,3, 6, 9, 15, 27 and here we only provide a brief introduction. The schematic of adaptive partitioning is illustrated in Fig. 2. The PAP potential is expressed in a many-body expansion formulation, where the energies of various QM/MM partitions are weighted and combined as follows:3 A *5, * A A A A A 𝑉 = 𝑉 A + ∑* &+, 𝑃& '𝑉& − 𝑉 ) + ∑&+, ∑-+&4, 𝑃& 𝑃- .𝑉&,- − '𝑉 + ∑1+&,-(𝑉1 − 𝑉 ))3 + *5, * A A A A A ∑*5: &+, ∑-+&4, ∑6+-4, 𝑃& 𝑃- 𝑃6 7𝑉&,-,6 − .𝑉 + ∑1+&,-,6 (𝑉1 − 𝑉 ) + ∑(8,9)+(&,-),(&,6),(-,6) .𝑉8,9 −

'𝑉 A + ∑1+8,9(𝑉1A − 𝑉 A ))33; + ⋯

(1)

Here, 𝑉 A is the QM/MM energy of the partitioning configuration (also called partition) with no buffer A group treated at the QM level, 𝑉&A with the i-th buffer group treated at the QM level, 𝑉&,with the i-th A and j-th buffer groups treated at the QM level, … 𝑉,,:,⋯,* with all N buffer groups at the QM level, and

𝑃& is the smoothing function for the i-th buffer group. Note that N can vary from one time step to another. We have chosen the smoothing function 𝑃& to be a fifth-order spline:2 𝑃(𝛼& ) = −6𝛼& ? + 15𝛼& B − 10𝛼& D + 1

(2)

where 𝛼& is a reduced distance for the i-th buffer group, 𝛼& = 1

1E 51min

(3)

FGH 51FIJ


7


Page 8 of 43

where distance 𝑟& = |𝐫& − 𝐫A | is between the buffer group at position 𝐫& and the active-zone center at position 𝐫A , and 𝑟max ≥ 𝑟𝑖 ≥ 𝑟min . The distance 𝑟& is group-based and is usually measured between representative atoms or the centers of mass of the groups. The smoothing function 𝑃& is continuous and differentiable over the domain [0, 1], and corresponds to all possible positions that a buffer group can take in the buffer zone.

Figure 2. Schematic of adaptive-partitioning QM/MM. The QM zone is centered at a selected molecule or ion A. The distance 𝑟& between a buffer group i and the QM-zone center satisfies 𝑟max ≥ 𝑟𝑖 ≥ 𝑟min . The QM, buffer, and MM groups are colored in orange, green, and blue, respectively.

While the PAP scheme is a Hamiltonian algorithm and conserves energy,3 nonphysical forces called “transition forces” arise from the derivative of the smoothing function 𝑃& .4, 9 The transition forces, if not negligible, may lead to erroneous solvation structures, e.g. by artificially forcing buffer groups out of the buffer zone.17 Various methods have been proposed to deal with the transition forces.4, 5, 9, 17, 27 In the mPAP scheme,9 external forces are introduced to exactly cancel out the transition forces. Conceptually, this means that the chemical potentials are equalized between different partitions, and consequently, the system described by mPAP is not a Hamiltonian one.21 The mPAP scheme conserves momentum and has been demonstrated to yield excellent structures and dynamics for systems simulated under constant temperatures through coupling to thermostats.9,

15, 21

Due to its simplicity and efficiency, we used the

mPAP9 scheme in this study.


8

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In this work, electrostatic embedding37, 133 is used for the QM/MM calculations, where the MM atomic charges of a large number of atoms surrounding the QM subsystem enter the effective QM Hamiltonians as one-electron operators. Although it is debatable how well the MM charges work as parameters for electrostatic embedding, this treatment offers the advantage that the QM wavefunction can be more realistic than that in mechanical embedding because it is polarized by the surrounding environment. Usually it is not practical or not necessary to enroll all MM atoms as background point charges in electrostatic embedding. Often, the background charges are selected using distance-based criteria, e.g. within a cut-off distance from a specific atom or spatial location. In our implementation of electrostaticembedding adaptive QM/MM, to avoid discontinuity in potentials, the background point charges are pre-selected and do not change except when treating buffer groups.6 Specifically, a group in the buffer zone will be treated as part of the background point charges in the QM/MM partitions where the group is not classified as QM. If a group exits the buffer zone and enters the environmental zone, it will always be treated as part of the background point charges. Currently, no other updates are allowed to the list of background point charges. 2.2. Extension of Proton Indicator Definition to Monatomic Ion and Carboxyl Group While the active-zone center can be conveniently tagged at an ordinary ion or molecule, here it is set to the H+ indicator.15 The position of the H+ indicator approximates the location of the excess H+. Developed and tested in our earlier works, the H+ indicator has been shown to be a faithful representation of the delocalized H+ in bulk water, yielding comparable results with other H+ representations such as the center of excess charge in dynamics simulations.15, 21, 57, 134 In this work, we generalized the H+ indicator definition because Cl– ion and carboxyl group are part of the H+ translocation path.


9


Page 10 of 43

Figure 3. Definition of H+ indicator exemplified by H+ exchange between a water molecule and a Cl– ion (green sphere), where D denotes the donor, HR the m-th H bonded to D, A- the j-th potential acceptor, and I the indicator (purple sphere). The red and silver spheres represent O and H atoms, respectively. The vector from D to A- is 𝐫DAT , and that from D to HR is 𝐫UVW . The enlisting distance 𝑟XYZ[ is indicated by the black dashed line and the blue circle centered at D.

The position of the H+ indicator is a linear combination of the coordinates of the donor D and all possible acceptors A- within a predefined enlisting radius 𝑟XYZ[ from D (Fig. 3): ,

𝐗 Y = ^ .𝐗 U + ∑d-+, ∑c R+, 𝑔'𝜌R- )𝐗 bT 3 _

(4)

Here, 𝐗 U and 𝐗 bT are the Cartesian coordinates of the donor oxygen and the j-th possible acceptor oxygen respectively, 𝑔'𝜌R- ) is the weight function associated with 𝜌R- (the ratio of the projected donor-acceptor vector), 𝑀 is the total number of HR (hydrogen atoms bonded to D), 𝐽 is the total number of possible acceptors within a radius of size 𝑟XYZ[ , and 𝑔Y is a normalization factor. The ratio 𝜌Ris a metric for how close a given atom HR is to donor D versus a possible acceptor A- and is calculated according to the following equation:15, 134, 135 𝜌R- =

𝐫ghW ∙ 𝐫gjT

(5)

l

k𝐫gj k T


10

Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where 𝐫UVW = 𝐗 VW −𝐗 U and 𝐫UbT = 𝐗 bT −𝐗 U . The linear combination of coordinates is normalized by 𝑔Y , which is calculated according to equation (6) using the weight function described in equation (7): 𝑔Y = 1 + ∑d-+, ∑c R+, 𝑔'𝜌R- )

(6)

0 if 1 ≤ 𝑥 ? B D 𝑔(𝑥) = n−6𝑥 + 15𝑥 − 10𝑥 + 1 if 0 ≤ 𝑥 < 1 1 if 𝑥 < 0

(7)

The weight function 𝑔'𝜌R- ) depends on the reduced variable x, which is now determined by t sWT 5sWT

𝑥 = 𝑥'𝜌R- ) = 1 − s

(8)

t max 5sWT

𝜌max = 1 t

t 1DH

(9)

t DH 41AH

u u u where 𝜌R= (𝑟UV /𝑟XYZ[ ) , 𝑟UV is a parameter that is set to slightly larger than the equilibrium D–H

bond distance (e.g. 1.00 Å in a hydronium ion) to reduce the sensitivity to the D–H vibrations near u equilibrium, 𝑟AH is similarly a parameter for the A–H bond distance, 𝑟XYZ[ is the threshold distance for

enlisting possible acceptors, and 𝜌max represents the percentage of the distance that a transferring H+ needs to travel from the donor towards an acceptor before the donor and acceptor switch status. In this u u work, the parameters 𝑟DH or 𝑟AH are based on gas-phase geometries optimized at the B3LYP/6-31G*

level for the protonated species (Table 1) except that for water the literature value15 is used. Note that in the original definition,15 𝜌max was simply fixed to 0.5 because only H+ relays between water molecules u u were considered, for which 𝑟DH = 𝑟AH . The purpose of using the enlisting threshold 𝑟XYZ[ (set to 3.5 Å)

is to gradually reduce the weight of possible acceptors as they move farther from the donor. As such, the indicator is dominated by the position of the donor and the most likely acceptor A, which is typically the acceptor closest to the donor.


11


Page 12 of 43

u u Table 1. Parameters 𝑟DH (or 𝑟AH ) for the H+ indicator in this work.

Donor or Acceptor

CHARMM Atom Type

u 𝑟DH [Å]

Water Oxygen

OT

1.0

Chloride Anion

CLA

1.4

Carboxylic Oxygen

OC

1.0

Tyrosine Oxygen

OH1

1.0

A donor-acceptor switch should occur at 𝜌R- = 𝜌max . However, because the simulations are conducted with finite numerical precision and finite time steps, the switch may actually occur at 𝜌R- > 𝜌max . Upon the D–A switch, the topology of the system is updated accordingly so that the involved MM calculations can be successfully carried out. To avoid discontinuities in the potential energy surface,37 the atom types and associated parameters (e.g. the MM atomic charges) do not change during the on-the-fly revision of molecular topology; only the topology is modified. As in our previous practice,15 the MM charges of all H atoms that are part of the H+ wire are set to the same as the MM charge of the H in a MM water molecule. This leads to a small, constant difference in the total charges of the QM subsystems in the associated MM and QM calculations, e.g. the QM subsystem may have a total net charge of +0.4 e in MM calculations but +1.0 e in QM calculations. We have so far not found any significant effects due to this difference. The revised description of H+ indicator can in principle describe tautomerization (H+ exchange between the two oxygen atoms) of the carboxyl group, although we did not observe such re-arrangement in the simulations in the work. 2.3. “Sushi-Roll” Partitions for One-Dimensional System In our previous adaptive QM/MM studies, an individual molecule or functional group comprised a group in adaptive partitioning. Such a strategy works well for highly diffusive systems, e.g. a solute in


12

Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bulk solvent, where any solvent molecule may independently diffuse into (or out of) the buffer zone or active zone. However, for transport through biological channels, the solute molecules being translocated along the pores are confined to movements approximately in one-dimension, and treating an individual molecule or functional group as a group may not be the optimal choice. First, the pore residues stay in about the same locations along the translocation path, even accounting for side chain flipping. Second, if an ordinary ion (e.g. Na+) or molecule is transported, water molecules in the pore often flow through together with the ion or molecule. On the other hand, if H+ is transported, it is probably the other extreme, where the water molecules serving as part of the H+ wire do not migrate as the excess H+ hops through as a charge defect or structural defect. In either case, the solvation shells of the solute are much more predictable, and this can be strategically exploited when defining the adaptive partitioning groups. Here we employ a “sushi-roll” strategy for ion transfer, where we slice the channel along the translocation path into individual groups (Fig. 4). Each group consists of multiple residues, ions, and/or water molecules. The pore model thus looks like a stack of sushi-rolls, which leads to the model’s name sushi-roll here, and readers should not confuse this with the ROLL136-138 algorithm, which was designed for anisotropic volume fluctuations and holonomic constraints. The dynamically relocated QM/MM boundary is treated by the redistributed charge scheme if it cuts through covalent bonds.6,

139

The

thickness and representative atom of each sushi-roll group should be tuned according to the desired sizes of the active and buffer zones (rmin and rmax) such that the active-zone is relayed smoothly along the channel while the number of buffer groups is limited to 1–2 at any given time step during the simulation. Because of the small number of buffer groups, the full PAP expansion often only contains terms up to the first or second order, thus reducing computational costs. The size of the active zone must be large enough so that the relevant chemical process can be computed with reasonably good accuracy, but not too large such that the calculations become unnecessarily expensive. For H+ transfer, our experience is that, a QM zone of ~5 Å radius will likely to yield good descriptions,21 and making each sushi-roll slice


13


Page 14 of 43

of ~5 Å thickness is a reasonable choice to begin with. In the case of an ordinary solute, where the water molecules surrounding the solute often travel with the solute, this strategy is useful as well, but these water molecules should be grouped together with the solute.

Figure 4. Illustrating the “sushi-roll” strategy for adaptive QM/MM simulations of ion (including H+) transport through a channel, as indicated by the red arrows. The pore residues are shown as sticks and individually colored. The channel is divided into five slices (groups 1–5) for adaptive partitioning, each containing one or more residues. The distance 𝑟& between each group and the transferred ion (A) is measured using one representative atom or the center of the mass (red dot) of the group. The buffer zone is defined by 𝑟max ≥ 𝑟𝑖 ≥ 𝑟min . The QM, buffer, and MM zones are highlighted in orange, green, and blue, respectively.

3. Computational Details 3.1. Model Preparation In a previous MD study by Han et al.,93 water molecules were observed spontaneously entering the hydrophobic pore of EcCLC to form a water wire with conformations that could facilitate H+ transfer between the internal (E203) and external (E148) gating residues. These water wires lasted for up to a


14

Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


few picoseconds before withdrawing from the pore. The model system used by Han et al.93 (denoted the Han model) guided our model construction in this work. We started by embedding an EcCLC protein dimer (PDB: 1OTS67) in a POPE lipid bilayer of 120 Å by 120 Å followed by adding two slabs of TIP3P140 water molecules of 30 Å thickness on either side of the lipid bilayer. To neutralize the system and to achieve a ~0.1 M physiological ion concentration, we replaced 151 randomly selected water molecules with 68 Na+ and 83 Cl– ions. To prepare the H+ wire, 6 water molecules were manually added in the pore. After superimposing our model with the Han model, the coordinates of the protein, Cl–cen, the 6 added water molecules, and the water molecules near the pore entrances (within 15 Å of these 6 water molecules) were set to those in the Han model. The final model system contains ~175,000 atoms. We then manually protonated the internal gating residue E203 in our model for two reasons. First, this is considered as one legitimate scenario for H+ transfer.109, 127 Second, by protonating E203, we created an downhill energy landscape for H+ migration all the way to E148.104 As a result, the H+ travels rapidly along the path, requiring quick tracking of the moving H+ and timely update of the molecular topology, QM/MM boundary, etc. in the adaptive partitioning. This presents a challenging situation for the testing of our method and code. After construction, the model system was minimized for 5,000 steps and subsequently equilibrated at the MM level with a time step of 1 fs for two subsequent runs for 5 ns and 10 ns under the NPT and NVT ensembles, respectively. All MM calculations in this work were performed using NAMD141 version 2.10. For the constant-pressure control, the Nosé-Hoover Langevin piston barostat142, 143 was used with a piston target set to 1 atm, period to 175 fs, and decay to 150 fs. For constant-temperature control, a Langevin thermostat144 was employed, with temperature set to 298 K and a damping coefficient set to 1 ps–1. After the first equilibration, the system reached its final dimensions of 119 Å × 119 Å × 118 Å in the x, y, and z dimensions respectively. During the second equilibration, we observed that one water


15


Page 16 of 43

molecule entered the pore from the intracellular bulk solution, leading to a total of 7 water molecules in the pore, which we call pore water molecules. To prepare for the H+ transfer, an excess H+ was added to the equilibrated model system by manually protonating the pore water molecule that is closest to the intracellular bulk solution. An additional short equilibration of 0.5 ns at the MM level was carried out to relax the hydronium ion, pore water molecules, and pore residues. All MM calculations used the CHARMM36145 force field for the protein, membrane, and ions, and the TIP3P water model. The parameters for the O–H bond, H–O–H angle, O and H charges, and O van der Waals interactions for the H3O+ were set to be the same as those for the TIP3P140 water. Periodic boundaries were described with the minimum image convention in order to be consistent with the later QM/MM calculations. For nonbonded interactions, a cutoff of 14 Å was utilized, with tapering switched at 13 Å. The pore water molecules, carboxyl O atoms of E203 and E148, and Cl–cen were restrained by harmonic potentials with force constants of 20 kcal×mol–1Å–2 to prevent significant deviations from the geometry in the Han model. These restraints were removed after the MM equilibration.


16

Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. (A) In conventional QM/MM simulations, the QM subsystem (shown as balls and thick sticks) consists of the Cl–cen ion, part of the E148 side chain, and 7 pore water molecules (W1–W7, positions taken from Ref. 93). The atoms surrounding the QM subsystem (thin sticks) are treated as background point charges in the electrostatic embedding. The lowest pore water W1 is manually protonated, and the indicator is displayed as a purple sphere (superimposed with the O atom of W1). Color code: C, cyan; H, white; O, red; N, blue; Cl, green. (B) Three sushi-roll models S1–S3 in adaptive QM/MM simulations. Only the QM atoms are shown, where the division of groups a, b, and c are indicated by dashed lines. The stars label the representative atoms of each group: the O atom of W1 for group a, of W5 for group b, and of W6 (model S1) or W7 (models S2 and S3) for group c.

3.2. Conventional QM/MM Dynamics Simulations As illustrated in Fig. 5A, the QM subsystem consisted of E148 side chain, Cl–cen, and the 7 pore water molecules, while the following pore residues are represented by background point charges in the electrostatic-embedding QM Hamiltonian: S107, I109, G146, R147, E148, G149, T151, V152, L186, A189, F190, F199, E203, F317, G355, I356, F357, A358, P359, L398, I402, A404, T407, Y445. This yielded 29 QM atoms and 321 MM atoms as background point charges. The H-link atom with the


17


redistributed charge scheme6,

139

Page 18 of 43

was employed to treat the QM/MM boundary passing through the

covalent bond connecting Cβ and Cγ of E148. The link atom (HL), which was assigned the HA2 atom type, was placed along the Cγ–Cβ bond with the bond length ratio

y'Cγ 5HL) y'Cγ 5Cβ )

= 0.72.139 The MM charge of

the Cβ atom was evenly redistributed to the mid-point of the Cβ–Ca and two Cβ–Hβ bonds. These redistributed point charges were added to the background point charge list. For all QM calculations in this work, the B3LYP146-148 functional was selected, with 6-31+G(d) basis sets chosen for Cl– and 631G(d) for all other atoms.149-153 Based on the subtractive definition of the QM/MM energy,37 the MM van der Waals interactions between a pair of atoms that both have a QM description in a QM/MM partition calculation should be exactly cancelled out in the MM calculations of the entire system and of the QM subsystem. However due to the finite precision in numerical calculations, there may be inexact cancellation between pairs of atoms that get very close to each other during a simulation. To avoid the inexact cancellations, we have set these van der Waals interactions to zero in our implementations by assigning new MM atom types to the QM atoms (e.g. OTQ for the water oxygen). In these new atom types, all parameters are the same as the corresponding MM atom types (e.g. OT for the water oxygen) except that the van der Waals parameters between QM atom pairs (OTQ–OTQ) are zeroed out using the NBFix term in the CHARMM parameter file. A Nosé-Hoover thermostat154 was employed with a coupling constant of 4 fs to retain the system at 298 K. A time step of 0.5 fs was utilized, and 40 trajectories were propagated from the same initial geometry but different atomic velocities. Each trajectory was propagated for up to 2 ps and was terminated early if completion of H+ transfer was observed. A smaller time step of 0.25 fs was also tested, but no significant difference was found in the trajectories. The trajectories were recorded at every 2.5 fs, and the H+ indicator position was calculated at every step. All QM/MM calculations in this work, including both conventional and adaptive ones, were performed with a local version of the QMMM155 code, which called NAMD141 version 2.10 for the


18

Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


involved MM calculations and Gaussian09156 for the QM calculations. The QMMM code synthesized the energies and gradients from both QM and MM calculations and propagated the trajectory accordingly. The simulations were carried out in the Comet supercomputer at the San Diego Supercomputer Center, and each trajectory was run on one Intel Haswell standard compute node consisting of two Intel Xeon E5-2680v3 processors at 2.5 GHz for 24 total cores, 128 GB of DRAM, and 320 GB SSD drive. 3.3. Adaptive QM/MM Dynamics Simulations The adaptive QM/MM simulations parameters were similar to the conventional ones except for the following differences. The QM subsystem in the conventional QM/MM simulations was divided into 3 sushi-roll groups. The buffer zone parameters 𝑟•‚ƒ and 𝑟•„… were set to 5.0 Å and 6.5 Å respectively. We have tested 3 different ways of dividing the pore into sushi rolls, which are denoted models S1 to S3 (Fig. 5B). As described in the methodology section, if the atoms in a group were not classified as QM atoms in an adaptive partition, they were included as background point charges in the embedded-QM calculations. The full PAP potential without truncation was employed, but as there were at most three (one active and two buffer) groups at any given time step, it was only possible to observe terms up to the 3rd order.


19


Page 20 of 43

4. Results and Discussion 4.1. Increase in Adaptive QM/MM Calculation Efficiency The sushi-roll method drastically increased the efficiency of adaptive calculations. During exploratory calculations where the one-molecule-per-AP-group scheme was used, the average wall-clock time per step was approximately 243 seconds. This is much slower than the sushi-roll scheme, where the average time per step was ~54 seconds. Although slower than the conventional QM/MM time per step of ~33 seconds, the sushi-roll time per step is expected to increase insignificantly for an even longer channel, because (i) the number of adaptive partitions and the maximum size of the QM subsystem in all these partitions will likely remain the same (containing no more than three sushi-roll groups) at any time step regardless of the total length of the channel and (ii) the time increases due to more background point charges in the embedded-QM calculations are modest. In contrast, the time per step needed in conventional QM/MM would grow dramatically, because the entire channel would have to be included as QM. 4.2. Flexibility of the H+ Transfer Pathway Although not as fluidic as bulk water molecules, the pore water molecules exhibited certain flexibility in our simulations, leading to multiple H+ transfer pathways. These pathways, which are tabulated in Table 2, can be characterized by the changes in the donors along the way. Here, we have excluded the back-and-forth fast H+ exchange within a donor−acceptor pair to yield the “net” transfer pathway for the system. For example, a transfer such as W1–W2–W1–W3–W4 will be recorded as W1–W3–W4 without W2. Although the numbers of the simulated trajectories are modest (40 for each tested QM/MM scheme), preventing us from making any bold claims, an interesting pattern emerged. In the conventional QM/MM calculations, one specific route dominated: Pathway 1, which is the most “direct” pathway W1–W3–W5–W6–W7–E148, accounted for about three quarters of the H+ transfer pathways. Two


20

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


routes tied for second place for the conventional calculations, each with 8% of the pathways: one where Cl–cen was included as a donor in the transport pathway and one where it was not. In our simulations, the Cl– ion became a donor if 𝑅(H′ − Cl) ≤ 1.4 Å, where H′ denotes the nearest H atom to Cl–cen at any time (the identity of H′ could change during a simulation). Pathways 1, 2, and 3 are illustrated in Fig. 6. It is encouraging to see that the adaptive QM/MM calculations reproduced the reference conventional QM/MM pathways quite well. All three sushi-roll models predicted Pathway 1 to be the most prominent route, representing 50–78% of the pathways. The agreement even extended to the second most popular route, Pathway 2, which accounted for ~8–13% of the total pathways. Overall, the agreement between the conventional and adaptive QM/MM simulations is, while not perfect, quite satisfactory.


21


Page 22 of 43

Table 2. Observed pathways for H+ transfer in conventional and adaptive QM/MM simulations. a ID Pathway

Conventional S1 S2 S3

1

W1—W3—W5—W6—W7—E148

29

20 29 31

2

W1—W3—W4—W5—W6—W7—E148

3

5

5

3

3

W1—W3—W5—Cl–cen—W6—W7—E148

3

0

1

1

4

W1—W3—W4—W6—W7—E148

2

4

0

2

5

W1—W3—W4—E148

1

0

0

0

6

W1—W3—W5—W4—W7—E148

1

0

0

0

7

W1—W3—W6—W7—E148

1

2

2

0

8

W1—W2—W4—W5—W6—W7—E148

0

1

0

2

9

W1—W3—W5—W4—W6—W7—E148

0

0

0

1

10 W1—W3—W5—W6—E148

0

3

1

0

11 W1—W3—W5—W7—E148

0

1

2

0

12 W1—W2—W4—W5—W6—E148

0

1

0

0

13 W1—W3—W5—W2—W6—E148

0

1

0

0

14 W1—W3—W4—W6—E148

0

1

0

0

15 W1—W2—W4—W6—W7—E148

0

1

0

0

40

40 40 40

Total a

S1, S2, and S3 denote the three sushi-roll models. See Fig. 5 for the labels W1 to W7 of pore water

molecules.


22

Page 23 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Scheme of the three most popular pathways for H+ transfer observed in conventional QM/MM simulations. Only the QM atoms are shown as balls and sticks: red, O; green, Cl; cyan, C; white, H.

4.3. A Two-Stage Journey of H+ Transfer to E148 Despite the flexibility in the H+ transfer pathways, visual inspection of the trajectories found that the overall pictures of H+ migration are qualitatively similar in all simulations. The general features are exemplified in Fig. 7 by the snapshots extracted from a representative trajectory of the conventional QM/MM simulations. First, we found that in general the excess H+ moved quickly (within 2 ps) to E148 and protonated it, in line with the overall downhill landscape in the potential surface of the scenario that we are simulating. Second, the H+ relay process could be largely divided into two stages. In the first stage, the H+ moved across the lower section of the H+ wire towards the Cl– ion and was then shared by the anion and its nearest water molecule (snapshot at 170 fs in Fig. 7). The H+ wandered near the Cl– ion for some time while the water molecules in the upper H+ wire aligned for transfer. The alignment included the bending of E148 towards the H+ and the changing of the positions and orientations of the nearby pore water molecules to create a hydrogen-bonding chain suitable for H+ relay between E203 and the Cl–. Once the path was in place, the H+ translocation entered the second stage, in which the H+ departed from the Cl– ion’s solvation shell and hopped rapidly (10 Å) distance crossed in the first stage needed more bridging water molecules. Due to the mobility of the water molecules, it is conceivably much less probable to have a path for nearly-concerted H+ transfer across a longer distance, because more water molecules will be required to be in specific positions and


25


Page 26 of 43

orientations at the same time. Therefore, H+ transport across longer distances is more likely to be dominated by step-wise movements and are slower. Therefore, it seems that the Cl–cen ion not only provides electrostatic attractions to help H+ migration toward the extracellular side but also eases the travel by serving as a “transit hub” that cuts a long journey into two successive shorter trips. Pre-organizing the H+ transfer pathway into two sections effectively reduces the conformational space of the water wire that the system needs to explore in order to move the H+ forward. This implies that the H+ transport will be faster via the two-stage transport assisted by Cl–cen than via a single-stage transport without such an assist, even though the mobility of the pore water molecules has already been significantly reduced compared with bulk water molecules.


26

Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 9. (A) Average velocity of the H+ indicator (I) advancing along the reaction coordinate 𝑅(I − O‘ ), i.e. the distance between the I and E148 Oε. S1–S3 denote the three adaptive-partitioning sushi-roll models. The velocity was computed via finite differentiation

Œy(Y5’“ ) Œ•

≈

y• (I5Oε )5y•–— (I5Oε ) ∆•

, where n indexes the snapshots, and ∆𝑡 is

the elapsed time between two successive snapshots. The velocities were binned in the range of 1 ≤ 𝑅(I − Oε ) ≤ 15 Å with a bin width of 0.5 Å according to the average value of the reaction coordinate over two successive snapshots

y•(I5Oε)4y•–— (I5Oε ) :

. The error bars indicate the standard deviation of each bin. (B) Average distance

between I and Cl–cen, 𝑅(I − Cl), as a function of 𝑅(I − O‘ ). The binning procedure is similar to that for

Œy(Y5’“ ) Œ•

,

and the error bars indicate the standard deviation of each bin.

4.4. Interactions between the H+ and Cl–cen One interesting question in the study of EcCLC is how the migrating H+ interacts with (or maybe even protonates) the bound Cl– ions. To get an insight into this problem, we have parametrically plotted the


27


Page 28 of 43

distance between the H+ indicator and Cl–cen, 𝑅(I − Cl) and the reaction coordinate 𝑅(I − O‘ ) as functions of simulation time in Fig. 10. The average curves, each over 40 trajectories of a given QM/MM scheme, are shown in Fig. 9B. These “V”-shape average curves are almost superimposable, implying that the H+…Cl– interactions are similar in the conventional and adaptive treatments. Starting from the upper-right corner, 𝑅(I − Cl) first decreases as 𝑅(I − O‘ ) shortens, until 𝑅(I − O‘ ) reaches 6–7 Å. Afterwards, 𝑅(I − Cl) increases as 𝑅(I − O‘ ) continues to reduce. While the average curves are smooth and show the same general trend, the individual curves (Fig. 10) are more rugged and display recrossings, reflecting the complexity of the H+ motion. The re-crossings are especially noticeable near the dip of the “V”, where the H+ stayed in the vicinity of the Cl– ion. Some individual curves even reach 𝑅(I − Cl) = 0, signifying that the Cl– ion had become the donor during the simulations. Indeed, we found that the Cl– ion became donor in ~53% of the conventional QM/MM trajectories and in 60%, 60%, and 68% of the S1, S2, and S3 adaptive QM/MM trajectories, respectively. Note that these percentages are calculated based on the H+ indicator position that was recorded every time step, so they are somewhat higher than what one will assume from Fig. 10, which is plotted based on the trajectories saved less frequently (every 5 steps). When serving as the donor during the H+ migration, Cl–cen participated mostly in the back-and-forth H+ exchange with a given pore water molecule and occasionally in the relay of the H+ to another water molecule.


28

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 10. Distance between the H+ indicator (I) and Cl–cen as a function of the distance between I and E148 Oε. S1– S3 denote the adaptive sushi-roll models S1–S3, respectively.

These data, together with our earlier analysis, indicated that the migrating H+ was shared by Cl–cen and nearby water molecules for a short but definite period of time before the second stages of the transport process. However, one may argue that the indicator is not a real atom. Therefore, we calculated the distance between Cl–cen and the nearest H atom (H′), 𝑅(H′ − Cl) for the trajectories and summarized the results in Table 3. Almost all trajectories had H atoms within 2 Å of Cl–cen, and more than 78% of the trajectories within 1.8 Å;157 the latter is the widely accepted radius of a Cl– ion. Of particular interest is the smallest 𝑅(H′ − Cl) of each trajectory. The average smallest 𝑅(H′ − Cl) over the 40 trajectories of a given QM/MM scheme are in the range of 1.52–1.60 Å, with the minimum smallest 𝑅(H′ − Cl) at 1.21– 1.26 Å. Compared with the gas-phase HCl equilibrium bond length of 1.275 Å,158, 159 it appears that protonation of the bound Cl– ion in EcCLC is possible, albeit only transiently, during the transport of the H+. Such a possibility had been previously suggested by experiments70 and computations.115 Of course, one must bear in mind that our observation is based on the specific model system and level of theory


29


Page 30 of 43

employed in this work, and the picture may somewhat change if different models and/or levels of theory are utilized. Moreover, as mentioned earlier, it has been shown that H+ can be transported under many scenarios including those in the absence of bound Cl–.61, 127 Thus, the present study alone will not be able to address the 2.2:1 stoichiometric ratio between the Cl– influx and H+ efflux, which is likely an “average” result of many contributing scenarios.127

Table 3. Statistics of 𝑅(H′ − Cl), the distance between Cl–cen and the nearest H atom (H′) a Percentage of Trajectories 𝑅(H′ − Cl) ≤ 2.0 𝑅(H′ − Cl) ≤ 1.8

Smallest 𝑅(H′ − Cl) b Average ± Standard Deviation Minimum

Conventional

98%

78%

1.60 ± 0.17

1.23

S1

98%

93%

1.52 ± 0.20

1.25

S2

100%

93%

1.55 ± 0.15

1.26

S3

100%

95%

1.53 ± 0.15

1.21

a

Based on saved trajectories. S1–S3 denote the three adaptive sushi-roll models. All distance in Å.

b

The smallest value of 𝑅(H′ − Cl) recorded in a trajectory.

Interestingly, in the adaptive sushi-roll QM/MM calculations, the migrating H+ seemed to come slightly closer on average to the bound Cl– ion than in the conventional QM/MM simulations, as judged from how many trajectories possessed 𝑅(H′ − Cl) ≤ 1.8 Å and from the average smallest 𝑅(H′ − Cl), both shown in Table 3. The subtle differences reflect the different QM subsystems at each time step between the adaptive and conventional treatments and between different sushi-roll models. In a sushiroll model, the central group (b) is very small, containing only the Cl– ion and 3–4 water molecule. Such a small group may somewhat overestimate the interactions between the Cl– ion and the solvated H+ water molecule. Fortunately, this group is not the sole QM contributing group, because the lower and


30

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


upper groups also contribute through the higher-order terms of the PAP potential. In the calculations of these higher-order terms, the central group is combined with the other groups to produce larger QM subsystems. Additionally, neighboring groups also polarize the QM wavefunctions of the central group through their representation as background point charges in the embedded-QM calculations of the central group. The contributions from other groups help to achieve a more realistic description that is closer to that by the larger QM subsystem in conventional QM/MM. How well the agreement is between the conventional and adaptive QM/MM will depend on the extent of electron delocalization between the ion and its surroundings. This underlines the importance of properly forming sushi-roll groups in adaptive QM/MM simulations. A sushi-roll slice of 5–6 Å thickness covers typically the first and second solvation shells of a small ion and is therefore recommended.

4. Conclusions: In this work, we carried out adaptive QM/MM dynamics simulations of H+ transfer through the CLC Cl–/H+ antiporter from E. coli. This is the first adaptive QM/MM simulation of ion transport through a channel and one of a few to use a high-level QM (DFT) method. Two technical advancements were made: (i) We have extended the H+ indicator definition so that it can be used for the H+ relay involving Cl– ions and carboxyl groups as H+ donors/acceptors. (ii) We have proposed a new “sushi-roll” scheme for setting up buffer groups in adaptive partitions for ion transport along one-dimensional channels. Although we have only simulated one particular scenario of H+ migration, the adaptive QM/MM reproduced the suite of conventional QM calculations quite satisfactorily in all aspects that we have examined, revealing essentially the same picture of H+ translocation and asserting the validity of the AP schemes. The computational costs of the adaptive QM/MM simulations are comparable with those of the conventional QM/MM for this relatively short channel, and it is reasonable to anticipate that the adaptive simulations will increasingly become more competitive when longer channels are simulated.


31


Page 32 of 43

We expect that the adaptive QM/MM algorithms can be applied to study the transport process in many other one-dimensional systems. Of course, adaptive QM/MM is not a replacement of conventional QM/MM. In many cases such as in umbrella sampling, it can be more convenient to divide the proton transport pathways into segments, which can be simulated separately by conventional QM/MM. Our simulations reveal a two-stage journey of the H+ migration in the presence of the Cl– ion bound at the central binding site (Cl–cen) of EcCLC. In the first stage, the H+ hops through the lower half of the pore to be near Cl–cen and stays there waiting for the water wire in the upper half of the pore to be ready. As soon as the upper half of the water wire is in position, the H+ transport enters the second stage, where the H+ moves quickly through the upper section of water wire to protonate the external gating residue E148. Cutting a long path for H+ transfer into two successive shorter sections with the assist by the bound Cl– ion reduces the conformational space of the water wire that the system needs to explore and should ease the H+ transport. The H+ relay path, consisting of pore water molecules, Cl–cen, and E148, exhibits certain mobility, and multiple pathways have been observed. In some of the simulated trajectories, the H+ came to very close distances (< 1.3 Å) from the bound Cl– ion, implying transient protonation of the Cl– ion, which had been suggested by previous experiments and computations.

ASSOCIATED CONTENT The supporting information (PDF file) contains the zeros of energy of adaptive partition groups and the snapshots of representative trajectories from adaptive QM/MM simulations. This information is available free of charge via the Internet at http://pubs.acs.org


32

Page 33 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


AUTHOR INFORMATION Corresponding Author Hai Lin, [email protected] Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Funding Sources This work is supported by the NSF (CHE-1564349), Camille & Henry Dreyfus Foundation (TH-14028), Research Corporation for Advancement (25793), and NVIDIA Corporation. C.M.G. and B.O.A are supported by the University of Colorado Denver through the Undergraduate Research Opportunity Program and the MARC U-STAR program. This work used XSEDE under grant CHE-140070, supported by NSF grant number ACI-1053575.

ACKNOWLEDGMENT We are grateful to Prof. Emad Tajkhorshid for making the simulated geometries of the water wire in EcCLC available to us.

ABBREVIATIONS EcCLC, Escheria coli Chloride Ion Channel; MM, molecular-mechanics; QM, quantum-mechanics; QM/MM, combined molecular-mechanics/quantum-mechanics; PAP, permuted adaptive partitioning


33


Page 34 of 43

REFERENCES 1. Kerdcharoen, T.; Liedl, K. R.; Rode, B. M., A QM/MM simulation method applied to the solution of Li in liquid ammonia. Chem. Phys. 1996, 211, 313-323. 2. Kerdcharoen, T.; Morokuma, K., ONIOM-XS: an extension of the ONIOM method for molecular simulation in condensed phase. Chem. Phys. Lett. 2002, 355, 257-262. 3. Heyden, A.; Lin, H.; Truhlar, D. G., Adaptive partitioning in combined quantum mechanical and molecular mechanical calculations of potential energy functions for multiscale simulations. J. Phys. Chem. B 2007, 111, 2231-2241. 4. Bulo, R. E.; Ensing, B.; Sikkema, J.; Visscher, L., Toward a practical method for adaptive QM/MM simulations. J. Chem. Theory Comput. 2009, 5, 2212-2221. 5. Nielsen, S. O.; Bulo, R. E.; Moore, P. B.; Ensing, B., Recent progress in adaptive multiscale molecular dynamics simulations of soft matter. Phys. Chem. Chem. Phys. 2010, 12, 12401-12414. 6. Pezeshki, S.; Lin, H., Adaptive-partitioning redistributed charge and dipole schemes for QM/MM dynamics simulations: On-the-fly relocation of boundaries that pass through covalent bonds. J. Chem. Theory Comput. 2011, 7, 3625-3634. 7. Takenaka, N.; Kitamura, Y.; Koyano, Y.; Nagaoka, M., The number-adaptive multiscale QM/MM molecular dynamics simulation: Application to liquid water. Chem. Phys. Lett. 2012, 524, 5661. 8. Bulo, R. E.; Michel, C.; Fleurat-Lessard, P.; Sautet, P., Multiscale modeling of chemistry in water: Are we there yet? J. Chem. Theory Comput. 2013, 9, 5567-5577. 9. Pezeshki, S.; Davis, C.; Heyden, A.; Lin, H., Adaptive-partitioning QM/MM dynamics simulations: 3. Solvent molecules entering and leaving protein binding sites J. Chem. Theory Comput. 2014, 10, 4765-4776. 10. Pezeshki, S.; Lin, H., Recent developments in QM/MM methods towards open-boundary multiscale simulations. Mol. Simulat. 2015, 41, 168-189. 11. Waller, M. P.; Kumbhar, S.; Yang, J., A density-based adaptive quantum mechanical/molecular mechanical method. ChemPhysChem 2014, 15, 3218-3225. 12. Watanabe, H. C.; Kubař, T.; Elstner, M., Size-consistent multipartitioning QM/MM: A stable and efficient adaptive QM/MM method. J. Chem. Theory Comput. 2014, 10, 4242-4252. 13. Böckmann, M.; Doltsinis, N. L.; Marx, D., Adaptive switching of interaction potentials in the time domain: An extended lagrangian approach tailored to transmute force field to QM/MM simulations and back. J. Chem. Theory Comput. 2015, 11, 2429-2439. 14. Jiang, T.; Boereboom, J. M.; Michel, C.; Fleurat-Lessard, P.; Bulo, R. E. Proton transfer in aqueous solution: Exploring the boundaries of adaptive QM/MM. In Quantum Modeling of Complex Molecular Systems, Rivail, J.-L.; Ruiz-Lopez, M.; Assfeld, X., Eds.; Springer International Publishing: Cham, 2015, pp 51-91. 15. Pezeshki, S.; Lin, H., Adaptive-partitioning QM/MM for molecular dynamics simulations: 4. Proton hopping in bulk water. J. Chem. Theory Comput. 2015, 11, 2398-2411. 16. Pezeshki, S.; Lin, H. Recent developments in adaptive QM/MM. In Quantum Modeling of Complex Molecular Systems, Rivail, J.-L.; Ruiz-Lopez, M.; Assfeld, X., Eds.; Springer: New York, NY, 2015, pp 93-113. 17. Boereboom, J. M.; Potestio, R.; Donadio, D.; Bulo, R. E., Toward Hamiltonian adaptive QM/MM: Accurate solvent structures using many-body potentials. J. Chem. Theory Comput. 2016, 12, 3441-3448. +


34

Page 35 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


18. Duster, A.; Garza, C.; Lin, H. Adaptive partitioning QM/MM dynamics simulations for substrate uptake, product release, and solvent exchange. In Computational Approaches for Studying Enzyme Mechanism, Voth, G. A., Ed.; Elsevier: Cambridge, MA, 2016, pp 342-358 19. Zheng, M.; Waller, M. P., Adaptive quantum mechanics/molecular mechanics methods. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2016, 6, 369-385. 20. Dohm, S.; Spohr, E.; Korth, M., Developing adaptive QM/MM computer simulations for electrochemistry. J. Comput. Chem. 2017, 38, 51-58. 21. Duster, A. W.; Wang, C.-H.; Garza, C. M.; Miller, D. E.; Lin, H., Adaptive quantum/molecular mechanics: what have we learned, where are we, and where do we go from here? Wiley Interdiscip. Rev. Comput. Mol. Sci. 2017, 7, e1310-n/a. 22. Field, M. J., An algorithm for adaptive QC/MM simulations. J. Chem. Theory Comput. 2017, 13, 2342-2351. 23. Zheng, M.; Kuriappan, J. A.; Waller, M. P., Toward more efficient density-based adaptive QM/MM methods. Int. J. Quantum Chem. 2017, 117, 1-8. 24. Boereboom, J. M.; Fleurat-Lessard, P.; Bulo, R. E., Explicit solvation matters: Performance of QM/MM solvation models in nucleophilic addition. J. Chem. Theory Comput. 2018, 14, 1841-1852. 25. Hofer, T. S.; Hünenberger, P. H., Absolute proton hydration free energy, surface potential of water, and redox potential of the hydrogen electrode from first principles: QM/MM MD free-energy simulations of sodium and potassium hydration. J. Chem. Phys. 2018, 148, 222814/1-28. 26. Zheng, M.; Waller, M. P., Yoink: An interaction-based partitioning API. J. Comput. Chem. 2018, 39, 799-806. 27. Duster, A.; Wang, C.-H.; Lin, H., Adaptive QM/MM for molecular dynamics simulations: 5. On the energy-conserved permuted adaptive-partitioning schemes. Molecules 2018, 23, 2170/1-16. 28. Warshel, A.; Levitt, M., Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227-249. 29. Singh, U. C.; Kollmann, P. A., A combined ab initio quantum mechanical and molecular mechanical method for carrying out simulations on complex molecular systems: Applications to the CH Cl + Cl exchange reaction and gas phase protonation of polyethers. J. Comput. Chem. 1986, 7, 718730. 30. Field, M. J.; Bash, P. A.; Karplus, M., A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations. J. Comput. Chem. 1990, 11, 700-733. 31. Gao, J., Methods and applications of combined quantum mechanical and molecular mechanical potentials. Rev. Comput. Chem. 1996, 7, 119-185. 32. Monard, G.; Merz, K. M., Jr., Combined quantum mechanical/molecular mechanical methodologies applied to biomolecular systems. Acc. Chem. Res. 1999, 32, 904-911. 33. Hammes-Schiffer, S., Theoretical perspectives on proton-coupled electron transfer reactions. Acc. Chem. Res. 2000, 34, 273-281. 34. Sherwood, P. Hybrid quantum mechanics/molecular mechanics approaches. In Modern Methods and Algorithms of Quantum Chemistry, Grotendorst, J., Ed.; John von Neumann-Instituts: Jülich, 2000; Vol. 3, pp 285-305. 35. Gao, J.; Truhlar, D. G., Quantum mechanical methods for enzyme kinetics. Annu. Rev. Phys. Chem. 2002, 53, 467-505. 36. Riccardi, D.; Schaefer, P.; Yang, Y.; Yu, H.; Ghosh, N.; Prat-Resina, X.; König, P.; Li, G.; Xu, D.; Guo, H., et al., Development of effective quantum mechanical/molecular mechanical (QM/MM) methods for complex biological process J. Phys. Chem. B 2006, 110, 6458-6469. −

3


35


Page 36 of 43

37. Lin, H.; Truhlar, D. G., QM/MM: What have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2007, 117, 185-199. 38. Senn, H. M.; Thiel, W., QM/MM methods for biological systems. Top. Curr. Chem. 2007, 268, 173-290. 39. Hu, H.; Yang, W., Free energies of chemical reactions in solution and in enzymes with ab initio quantum mechanics/molecular mechanics methods. Annu. Rev. Phys. Chem. 2008, 59, 573-601. 40. Sherwood, P.; Brooks, B. R.; Sansom, M. S. P., Multiscale methods for macromolecular simulations. Curr. Opin. Struct. Biol. 2008, 18, 630-640. 41. Sabin, J. R.; Brändas, E., Combining quantum mechanics and molecular mechanics. Some recent progresses in QM/MM methods. Academic Press: San Diego, CA, 2010; p 1-416. 42. Ferrer, S.; Ruiz-Pernia, J.; Marti, S.; Moliner, V.; Tunon, I.; Bertran, J.; Andres, J. Hybrid schemes based on quantum mechanics/molecular mechanics simulations: Goals to success, problems, and perspectives. In Advances in Protein Chemistry and Structural Biology, Vol 85: Computational Chemistry Methods in Structural Biology, Christov, C., Ed.; Elsevier Academic Press Inc: San Diego, CA, 2011; Vol. 85, pp 81-142. 43. Menikarachchi, L. C.; Gascon, J. A., QM/MM approaches in medicinal chemistry research. Curr. Top. Med. Chem. 2010, 10, 46-54. 44. Wallrapp, F. H.; Guallar, V., Mixed quantum mechanics and molecular mechanics methods: Looking inside proteins. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 315-322. 45. Woodcock, H. L.; Miller, B. T.; Hodoscek, M.; Okur, A.; Larkin, J. D.; Ponder, J. W.; Brooks, B. R., MSCALE: A general utility for multiscale modeling. J. Chem. Theory Comput. 2011, 7, 12081219. 46. Chung, L. W.; Hirao, H.; Li, X.; Morokuma, K., The ONIOM method: its foundation and applications to metalloenzymes and photobiology. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 327-350. 47. Lonsdale, R.; Harvey, J. N.; Mulholland, A. J., A practical guide to modelling enzyme-catalysed reactions. Chem. Soc. Rev. 2012, 41, 3025-3038. 48. Monari, A.; Rivail, J.-L.; Assfeld, X., Theoretical modeling of large molecular systems. Advances in the local self consistent field method for mixed quantum mechanics/molecular mechanics calculations. Acc. Chem. Res. 2012, 46, 596-603. 49. Wu, R.; Cao, Z.; Zhang, Y., Computational simulations of zinc enzyme: Challenges and recent advances. Prog. Chem. 2012, 24, 1175-1184. 50. Mennucci, B., Modeling environment effects on spectroscopies through QM/classical models. Phys. Chem. Chem. Phys. 2013, 15, 6583-6594. 51. Meier, K.; Choutko, A.; Dolenc, J.; Eichenberger, A. P.; Riniker, S.; van Gunsteren, W. F., Multi-resolution simulation of biomolecular systems: A review of methodological issues. Angew. Chem. Int. Ed. 2013, 52, 2820-2834. 52. Duarte, F.; Amrein, B. A.; Blaha-Nelson, D.; Kamerlin, S. C. L., Recent advances in QM/MM free energy calculations using reference potentials. BBA-Gen. Subjects 2015, 1850, 954-965. 53. Agmon, N., The Grotthuss mechanism. Chem. Phys. Lett. 1995, 244, 456-462. 54. Marx, D., Proton transfer 200 years after von Grotthuss: Insights from ab initio simulations. ChemPhysChem 2006, 7, 1848-1870. 55. Swanson, J. M. J.; Maupin, C. M.; Chen, H.; Petersen, M. K.; Xu, J.; Wu, Y.; Voth, G. A., Proton solvation and transport in aqueous and biomolecular systems: Insights from computer simulations. J. Phys. Chem. B 2007, 111, 4300-4314.


36

Page 37 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


56. Knight, C.; Voth, G. A., The Curious Case of the Hydrated Proton. Acc. Chem. Res. 2012, 45, 101-109. 57. Duster, A. W.; Lin, H., Restrained proton indicator in combined quantum-mechanics/molecularmechanics dynamics simulations of proton transfer through a carbon nanotube. J. Phys. Chem. B 2017, 121, 8585-8592. 58. Riccardi, D.; Konig, P.; Prat-Resina, X.; Yu, H. B.; Elstner, M.; Frauenheim, T.; Cui, Q., "Proton holes" in long-range proton transfer reactions in solution and enzymes: A theoretical analysis. J. Am. Chem. Soc. 2006, 128, 16302-16311. 59. Konig, P. H.; Ghosh, N.; Hoffmann, M.; Elstner, M.; Tajkhorshid, E.; Frauenheim, T.; Cui, Q., Toward theoretical analyis of long-range proton transfer kinetics in biomolecular pumps. J. Phys. Chem. A 2006, 110, 548-563. 60. Liang, R.; Swanson, J. M. J.; Madsen, J. J.; Hong, M.; DeGrado, W. F.; Voth, G. A., Acid activation mechanism of the influenza A M2 proton channel. Proc. Natl. Acad. Sci. 2016, 113, E6955E6964. 61. Lee, S.; Swanson, J. M. J.; Voth, G. A., Multiscale simulations reveal key aspects of the proton transport mechanism in the ClC-ec1 antiporter. Biophys. J. 2016, 110, 1334-1345. 62. Day, T. J. F.; Soudackov, A. V.; Čuma, M.; Schmitt, U. W.; Voth, G. A., A second generation multistate empirical valence bond model for proton transport in aqueous systems. J. Chem. Phys. 2002, 117, 5839-5849. 63. Peng, Y.; Swanson, J. M. J.; Kang, S.-g.; Zhou, R.; Voth, G. A., Hydrated excess protons can create their own water wires. J. Phys. Chem. B 2015, 119, 9212-9218. 64. Liang, R.; Li, H.; Swanson, J. M. J.; Voth, G. A., Multiscale simulation reveals a multifaceted mechanism of proton permeation through the influenza A M2 proton channel. Proc. Natl. Acad. Sci. 2014, 111, 9396-9401. 65. Chakrabarti, N.; Tajkhorshid, E.; Roux, B.; Pomès, R., Molecular basis of proton blockage in aquaporins. Structure 2004, 12, 65-74. 66. Maduke, M.; Miller, C.; Mindell, J. A., A decade of ClC chloride channels: Structure, mechanism, and many unsettled questions. Annu. Rev. Biophys. Biomol. Struct. 2000, 29, 411-438. 67. Dutzler, R.; Campbell, E. B.; MacKinnon, R., Gating the selectivity filter in ClC chloride channels. Science 2003, 300, 108-112. 68. Accardi, A.; Miller, C., Secondary active transport mediated by a prokaryotic homologue of ClC Cl channels. Nature 2004, 427, 803-807. 69. Matulef, K.; Maduke, M., The CLC 'chloride channel' family: revelations from prokaryotes. Mol. Membr. Biol. 2007, 24, 342-350. 70. Miller, C.; Nguitragool, W., A provisional transport mechanism for a chloride channel-type Cl /H exchanger. Phil. Trans. R. Soc. B 2009, 364, 175-180. 71. Jentsch, T. J., CLC chloride channels and transporters: From genes to protein structure, pathology and physiology. Crit. Rev. Biochem. Mol. Biol. 2008, 43, 3-36. 72. Accardi, A., Structure and gating of CLC channels and exchangers. J. Physiol. 2015, 593, 41294138. 73. Steinmeyer, K.; Ortland, C.; Jentsch, T. J., Primary structure and functional expression of a developmentally regulated skeletal muscle chloride channel. Nature 1991, 354, 301. 74. Dutzler, R., The ClC family of chloride channels and transporters. Curr. Opin. Struct. Biol. 2006, 16, 439-446. 75. Accardi, A.; Picollo, A., CLC channels and transporters: Proteins with borderline personalities. BBA-Biomembranes 2010, 1798, 1457-1464. −

−

+


37


Page 38 of 43

76. Jentsch, T. J., Discovery of CLC transport proteins: cloning, structure, function and pathophysiology. J. Physiol. 2015, 593, 4091-4109. 77. Picollo, A.; Pusch, M., Chloride/proton antiporter activity of mammalian CLC proteins ClC-4 and ClC-5. Nature 2005, 436, 420-423. 78. Scheel, O.; Zdebik, A. A.; Lourdel, S.; Jentsch, T. J., Voltage-dependent electrogenic chloride/proton exchange by endosomal CLC proteins. Nature 2005, 436, 424-427. 79. Feng, L.; Campbell, E. B.; MacKinnon, R., Molecular mechanism of proton transport in CLC Cl /H exchange transporters. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 11699-11704. 80. Picollo, A.; Xu, Y.; Johner, N.; Bernèche, S.; Accardi, A., Synergistic substrate binding determines the stoichiometry of transport of a prokaryotic H /Cl exchanger. Nat. Struct. Mol. Biol. 2012, 19, 525-531. 81. Accardi, A.; Walden, M.; Nguitragool, W.; Jayaram, H.; Williams, C.; Miller, C., Separate ion pathways in a Cl /H exchanger. J. Gen. Physiol. 2005, 126, 563-570. 82. Dutzler, R.; Campbell, E. B.; Cadene, M.; Chait, B. T.; MacKinnon, R., X-ray structure of a ClC chloride channel at 3.0 Å reveals the molecular basis of anion selectivity. Nature 2002, 415, 287-294. 83. Accardi, A.; Lobet, S.; Williams, C.; Miller, C.; Dutzler, R., Synergism between halide binding and proton transport in a CLC-type exchanger. J. Mol. Biol. 2006, 362, 691-699. 84. Nguitragool, W.; Miller, C., Uncoupling of a CLC Cl /H exchange transporter by polyatomic anions. J. Mol. Biol. 2006, 362, 682-690. 85. Walden, M.; Accardi, A.; Wu, F.; Xu, C.; Williams, C.; Miller, C., Uncoupling and turnover in a Cl /H exchange transporter. J. Gen. Physiol. 2007, 129, 317-329. 86. Jayaram, H.; Accardi, A.; Wu, F.; Williams, C.; Miller, C., Ion permeation through a Cl selective channel designed from a CLC Cl /H exchanger. Proc. Natl. Acad. Sci. 2008, 105, 1119411199. 87. Lim, H.-H.; Miller, C., Intracellular proton-transfer mutants in a ClC Cl /H exchanger. J. Gen. Physiol. 2009, 133, 131-138. 88. Elvington, S. M.; Liu, C. W.; Maduke, M. C., Substrate-driven conformational changes in ClCec1 observed by fluorine NMR. EMBO J. 2009, 28, 3090-3102. 89. Picollo, A.; Malvezzi, M.; Houtman, J. C. D.; Accardi, A., Basis of substrate binding and conservation of selectivity in the CLC family of channels and transporters. Nat. Struct. Mol. Biol. 2009, 16, 1294-1301. 90. Robertson, J. L.; Kolmakova-Partensky, L.; Miller, C., Design, function and structure of a monomeric ClC transporter. Nature 2010, 468, 844-847. 91. Lim, H.-H.; Shane, T.; Miller, C., Intracellular proton access in a Cl /H antiporter PLoS Biol 2012, 10, e1001441. 92. Basilio, D.; Noack, K.; Picollo, A.; Accardi, A., Conformational changes required for H /Cl exchange mediated by a CLC transporter. Nat. Struct. Mol. Biol. 2014, 21, 456-464. 93. Han, W.; Cheng, R. C.; Maduke, M. C.; Tajkhorshid, E., Water access points and hydration pathways in CLC H /Cl transporters. Proc. Natl. Acad. Sci. 2014, 111, 1819-1824. 94. Khantwal, C. M.; Abraham, S. J.; Han, W.; Jiang, T.; Chavan, T. S.; Cheng, R. C.; Elvington, S. M.; Liu, C. W.; Mathews, I. I.; Stein, R. A., et al., Revealing an outward-facing open conformational state in a CLC Cl /H exchange transporter. eLife 2016, 5, e11189/1-30. 95. Park, E.; Campbell, E. B.; MacKinnon, R., Structure of a CLC chloride ion channel by cryoelectron microscopy. Nature 2017, 541, 500-505. 96. Vien, M.; Basilio, D.; Leisle, L.; Accardi, A., Probing the conformation of a conserved glutamic acid within the Cl pathway of a CLC H /Cl exchanger. J. Gen. Physiol. 2017, 149, 523-529. −

+

+

−

−

+

−

−

+

+

−

−

+

−

−

+

+

+

+

-

−

−

+

−

+

−


38

Page 39 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


97. Cohen, J.; Schulten, K., Mechanism of anionic conduction across ClC. Biophys. J. 2004, 86, 836845. 98. Corry, B.; O'Mara, M.; Chung, S.-H., Conduction mechanisms of chloride ions in ClC-type channels. Biophys. J. 2004, 86, 846-860. 99. Miloshevsky, G. V.; Jordan, P. C., Anion pathway and potential energy profiles along curvilinear bacterial ClC Cl pores: Electrostatic effects of charged residues. Biophys. J. 2004, 86, 825-835. 100. Yin, J.; Kuang, Z.; Mahankali, U.; Beck, T. L., Ion transit pathways and gating in ClC chloride channels. Proteins: Struct., Funct., Bioinf. 2004, 57, 414-421. 101. Faraldo-Gómez, J. D.; Roux, B., Electrostatics of ion stabilization in a ClC chloride channel homologue from Escherichia coli. J. Mol. Biol. 2004, 339, 981-1000. 102. Bostick, D. L.; Berkowitz, M. L., Exterior site occupancy infers chloride-induced proton gating in a prokaryotic homolog of the ClC chloride channel. Biophys. J. 2004, 87, 1686-1696. 103. Bisset, D.; Corry, B.; Chung, S.-H., The fast gating mechanism in ClC-0 channels. Biophys. J. 2005, 89, 179-186. 104. Gervasio, F. L.; Parrinello, M.; Ceccarelli, M.; Klein, M. L., Exploring the gating mechanism in the ClC chloride channel via metadynamics. J. Mol. Biol. 2006, 361, 390-398. 105. Cheng, M. H.; Mamonov, A. B.; Dukes, J. W.; Coalson, R. D., Modeling the fast gating mechanism in the ClC-0 chloride channel. J. Phys. Chem. B 2007, 111, 5956-5965. 106. Engh, A. M.; Faraldo-Gómez, J. D.; Maduke, M., The role of a conserved lysine in chloride- and voltage-dependent ClC-0 fast gating. J. Gen. Physiol. 2007, 130, 351-363. 107. Engh, A. M.; Faraldo-Gómez, J. D.; Maduke, M., The mechanism of fast-gate opening in ClC-0. J. Gen. Physiol. 2007, 130, 335-349. 108. Kuang, Z.; Mahankali, U.; Beck, T. L., Proton pathways and H /Cl stoichiometry in bacterial chloride transporters. Proteins: Struct., Funct., Bioinf. 2007, 68, 26-33. 109. Wang, D.; Voth, G. A., Proton transport pathway in the ClC Cl /H antiporter. Biophys. J. 2009, 97, 121-131. 110. Coalson, R. D.; Cheng, M. H., Discrete-state representation of ion permeation coupled to fast gating in a model of ClC chloride channels: Comparison to multi-ion continuous space brownian dynamics simulations. J. Phys. Chem. B 2010, 114, 1424-1433. 111. Ko, Y. J.; Jo, W. H., Chloride ion conduction without water coordination in the pore of ClC protein. J. Comput. Chem. 2010, 31, 603-611. 112. Ko, Y. J.; Jo, W. H., Secondary water pore formation for proton transport in a CIC exchanger revealed by an atomistic molecular-dynamics simulation. Biophys. J. 2010, 98, 2163-2169. 113. Miloshevsky, G. V.; Hassanein, A.; Jordan, P. C., Antiport mechanism for Cl /H in ClC-ec1 from normal-mode analysis. Biophys. J. 2010, 98, 999-1008. 114. Coalson, R. D.; Cheng, M. H., Discrete-state representation of ion permeation coupled to fast gating in a model of CLC-chloride channels: Analytic estimation of the state-to-state rate constants. J. Phys. Chem. A 2011, 115, 9633-9642. 115. Kieseritzky, G.; Knapp, E.-W., Charge transport in the ClC-type chloride-proton anti-porter from Escherichia coli. J. Biol. Chem. 2011, 286, 2976-2986. 116. Smith, M.; Lin, H., Charge delocalization upon chloride ion binding in ClC chloride ion channels/transporters. Chem. Phys. Lett. 2011, 502, 112-117. 117. Zhang, Y.; Voth, Gregory A., The coupled proton transport in the ClC-ec1 Cl /H antiporter. Biophys. J. 2011, 101, L47-L49. 118. Cheng, M. H.; Coalson, R. D., Molecular dynamics investigation of Cl and water transport through a eukaryotic CLC transporter. Biophys. J. 2012, 102, 1363-1371. −

+

−

−

+

−

+

−

+

−


39


Page 40 of 43

119. Yu, T.; Wang, X.-Q.; Sang, J.-P.; Pan, C.-X.; Zou, X.-W.; Chen, T.-Y.; Zou, X., Influences of mutations on the electrostatic binding free energies of chloride ions in Escherichia coli ClC. J. Phys. Chem. B 2012, 116, 6431-6438. 120. Church, J.; Pezeshki, S.; Davis, C.; Lin, H., Charge transfer and polarization for chloride ions bound in ClC transport proteins: Natural bond orbital and energy decomposition analyses. J. Phys. Chem. B 2013, 117, 16029-16043. 121. Nieto-Delgado, P. G.; Arreola, J.; Guirado-López, R. A., Atomic charges of Cl ions confined in a model Escherichia coli ClC−Cl /H ion exchanger: a density functional theory study. Mol. Phys. 2013, 111, 3218-3233. 122. Chen, Z.; Beck, T. L., Free energies of ion binding in the bacterial CLC-Ec1 chloride transporter with implications for the transport mechanism and selectivity. J. Phys. Chem. B 2016, 120, 3129-3139. 123. Jiang, T.; Han, W.; Maduke, M.; Tajkhorshid, E., Molecular basis for differential anion binding and proton coupling in the Cl /H exchanger CLC-ec1. J. Am. Chem. Soc. 2016, 138, 3066-3075. 124. Lee, S.; Liang, R.; Voth, G. A.; Swanson, J. M. J., Computationally efficient multiscale reactive molecular dynamics to describe amino acid deprotonation in proteins. J. Chem. Theory Comput. 2016, 12, 879-891. 125. Lee, S.; Mayes, H. B.; Swanson, J. M. J.; Voth, G. A., The origin of coupled chloride and proton transport in a Cl /H antiporter. J. Am. Chem. Soc. 2016, 138, 14923-14930. 126. Chenal, C.; Gunner, M. R., Two Cl ions and a glu compete for a helix cage in the CLC proton/Cl antiporter. Biophys. J. 2017, 113, 1025-1036. 127. Mayes, H. B.; Lee, S.; White, A. D.; Voth, G. A.; Swanson, J. M. J., Multiscale kinetic modeling reveals an ensemble of Cl /H exchange pathways in ClC-ec1 antiporter. J. Am. Chem. Soc. 2018, 140, 1793-1804. 128. Wang, C.-H.; Duster, A. W.; Aydintug, B. O.; Zarecki, M. G.; Lin, H., Chloride ion transport by the E. coli CLC Cl /H antiporter: A combined quantum-mechanical and molecular-mechanical study. Frontiers in Chemistry 2018, 6, 62/1-16. 129. Zwanzig, R., From classical dynamics to continuous time random walks. J Stat Phys 1983, 30, 255-262. 130. Noé, F., Probability distributions of molecular observables computed from Markov models. J. Chem. Phys. 2008, 128, 244103/1-13. 131. Wang, W.; Cao, S.; Zhu, L.; Huang, X., Constructing Markov State Models to elucidate the functional conformational changes of complex biomolecules. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2018, 8, e1343. 132. Husic, B. E.; Pande, V. S., Markov state models: From an art to a science. J. Am. Chem. Soc. 2018, 140, 2386-2396. 133. Bakowies, D.; Thiel, W., Hybrid models for combined quantum mechanical and molecular mechanical approaches. J. Phys. Chem. 1996, 100, 10580-10594. 134. Wu, X.; Thiel, W.; Pezeshki, S.; Lin, H., Specific reaction path Hamiltonian for proton transfer in water: Reparameterized semiempirical models. J. Chem. Theory Comput. 2013, 9, 2672-2686. 135. Hofer, T. S.; Hitzenberger, M.; Randolf, B. R., Combining a dissociative water model with a hybrid QM/MM approach-a simulation strategy for the study of proton transfer reactions in solution. J. Chem. Theory Comput. 2012, 8, 3586-3595. 136. Yu, T.-Q.; Alejandre, J.; López-Rendón, R.; Martyna, G. J.; Tuckerman, M. E., Measurepreserving integrators for molecular dynamics in the isothermal–isobaric ensemble derived from the Liouville operator. Chem. Phys. 2010, 370, 294-305. −

−

–

–

+

+

+

-

–

−

+

+


40

Page 41 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


137. Tuckerman, M. E.; Alejandre, J.; López-Rendón, R.; Jochim, A. L.; Martyna, G. J., A Liouvilleoperator derived measure-preserving integrator for molecular dynamics simulations in the isothermal– isobaric ensemble. J. Phys. A: Math. Gen. 2006, 39, 5629. 138. Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L., Explicit reversible integrators for extended systems dynamics. Mol. Phys. 1996, 87, 1117-1157. 139. Lin, H.; Truhlar, D. G., Redistributed charge and dipole schemes for combined quantum mechanical and molecular mechanical calculations. J. Phys. Chem. A 2005, 109, 3991-4004. 140. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926-935. 141. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K., Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781-1802. 142. Martyna, G. J.; Tobias, D. J.; Klein, M. L., Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177-4189. 143. Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R., Constant pressure molecular dynamics simulation: The Langevin piston method. J. Chem. Phys. 1995, 103, 4613-4621. 144. Brünger, A.; Brooks Iii, C. L.; Karplus, M., Stochastic boundary conditions for molecular dynamics simulations of ST2 water. Chem. Phys. Lett. 1984, 105, 495-500. 145. Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D., Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ϕ, ψ and side-chain χ1 and χ2 dihedral angles. J. Chem. Theory Comput. 2012, 8, 3257-3273. 146. Becke, A. D., Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098-3100. 147. Becke, A. D., Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648-5652. 148. Lee, C.; Yang, W.; Parr, R. G., Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter 1988, 37, 785-789. 149. Ditchfield, R.; Hehre, W. J.; Pople, J. A., Self-consistent molecular-orbital methods. IX. An extended Gaussian-type basis for molecular-orbital studies of organic molecules. J. Chem. Phys. 1971, 54, 724-728. 150. Hehre, W. J.; Ditchfield, R.; Pople, J. A., Self-consistent molecular orbital methods. XII. Further extensions of Gaussian-type basis sets for use in molecular orbital studies of organic molecules. J. Chem. Phys. 1972, 56, 2257-2261. 151. Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; DeFrees, D. J.; Pople, J. A.; Gordon, M. S., Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. J. Chem. Phys. 1982, 77, 3654-3665. 152. Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R., Efficient diffuse functionaugmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li-F. J. Comput. Chem. 1983, 4, 294-301. 153. Frisch, M. J.; Pople, J. A.; Binkley, J. S., Quadratic configuration interaction. A general technique for determining electron correlation energies. J. Chem. Phys. 1984, 80, 3265-3269. 154. Hoover, W. G., Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695-1697. 155. Lin, H.; Zhang, Y.; Pezeshki, S.; Duster, A.; Truhlar, D. G. QMMM, Version 2.2.8.CO; University of Minnesota: Minneapolis, 2017.


41


Page 42 of 43

156. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., et al. Guassian09, Version B.01; Gaussian, Inc.: Wallingford, CT 2010. 157. Shannon, R., Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides. Acta Crystallographica Section A 1976, 32, 751-767. 158. Bunker, P. R., On the breakdown of the Born-Oppenheimer approximation for a diatomic molecule. J. Mol. Spectrosc. 1972, 42, 478-494. 159. Watson, J. K. G., The isotope dependence of the equilibrium rotational constants in 1Σ states of diatomic molecules. J. Mol. Spectrosc. 1973, 45, 99-113.


42

Page 43 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic


43

MM for Molecular Dynamics Simulations: 6

Recommend Documents