What Can and Cannot Be Learned from Molecular Dynamics

Dec 13, 2016 - This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (...
0 downloads 7 Views 4MB Size
Article pubs.acs.org/JPCB

What Can and Cannot Be Learned from Molecular Dynamics Simulations of Bacterial Proton-Coupled Oligopeptide Transporter GkPOT? Kalyan Immadisetty, Jeevapani Hettige, and Mahmoud Moradi* Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, United States S Supporting Information *

ABSTRACT: We have performed an extensive set of all-atom molecular dynamics (MD) simulations of a bacterial proton-coupled oligopeptide transporter (POT) in an explicit membrane environment. We have characterized both the local and global conformational dynamics of the transporter upon the proton and/or substrate binding, within a statistical framework. Our results reveal a clearly distinct behavior for local conformational dynamics in the absence and presence of the proton at the putative proton binding residue E310. Particularly, we find that the substrate binding conformation is drastically different in the two conditions, where the substrate binds to the protein in a lateral/vertical manner, in the presence/absence of the proton. We do not observe any statistically significant distinctive behavior in terms of the global conformational changes in different simulation conditions, within the time scales of our simulations. Our extensive simulations and analyses call into question the implicit assumption of many MD studies that local conformational changes observed in short simulations could provide clues to the global conformational changes that occur on much longer time scales. The linear regression analysis of quantities associated with the global conformational fluctuations, however, provides an indication of a mechanism involving the concerted motion of the transmembrane helices, consistent with the rocker-switch mechanism.

1. INTRODUCTION Proton-coupled oligopeptide transporters (POTs) are membrane proteins that belong to the major facilitator superfamily (MFS) of transporters.1−3 The MFS transporters including uniporters, antiporters, and symporters have extensively been studied structurally over the past decade.4−17 POTs are among the symporter members of the MFS, which import small peptides and peptide-like molecules18 across the cell membrane (e.g., in human small intestine) using the inwardly directed proton flow (Figure 1). A key feature among POTs is substrate promiscuity,19 which is partly attributed to the ability of the binding site to accommodate peptides in multiple orientations (e.g., vertical and lateral).20 Biophysical studies of POT transporter YjdL21 suggest the preference of the transporter for dipeptides over tripeptides. Substrate specificity/promiscuity in POTs is of great interest, particularly due to its potential relevance to biomedical applications.21−23 Human POT transporters PepT1 and PepT2 play a key role in absorbing and retaining dietary proteins.24 They recognize several important families of peptide-like drug compounds such as βlactam antibiotics25 and can improve the uptake of poorly absorbed/retained medications if attached to amino acids or dipeptides.26−28 PepT1 and PepT2 can uptake these so-called prodrugs, which may provide an important tool in the development of more effective medications.24 Mammalian POTs are yet to be crystallized; however, several bacterial POT members have recently been crystallized including PepTso,8,20 PepTst,9 PepTso2,23 and GkPOT.29 © 2016 American Chemical Society

Figure 1. Schematic representation of putative transport cycle in POTs. It is believed that the global transition occurs in either the absence (left) or presence (right) of both the proton (H+) and the substrate, while proton binding stabilizes the protein in an either OF (top) or IF (bottom) state to allow for substrate binding/unbinding.

Special Issue: Klaus Schulten Memorial Issue Received: September 26, 2016 Revised: December 13, 2016 Published: December 13, 2016 3644

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B Particularly, the POT transporter from Geobacillus kaustophilus bacterium (GkPOT) has been crystallized in a lipidic environment to a remarkably high resolution of 1.9−2.4 Å, in both the apo and substrate-bound states, and for both the wildtype protein and E310Q mutant, the latter mimicking the protonation of one of the putative proton binding sites, i.e., the E310 residue.29 The study of bacterial POT transporters could shed light on the transport mechanism of human PepT1 and PepT2, due to the high degree of sequence conservation within the transmembrane (TM) helices in POTs.24 While crystal structures provide a unique opportunity for establishing a structural basis for the mechanism of POT transporters, many mechanistic details remain uncharacterized, partly due to the elusiveness of the outward-facing (OF) state in POTs. All POT crystal structures, reported to date, are in an inward-facing (IF) state.8,9,20,23 In order to actively transport materials across the membrane, however, membrane transporters have to undergo global conformational changes to alternate between distinct IF and OF states (Figure 1), during which the substrate binding site is open only to one side of the membrane at any given moment. The “alternating access” mechanism described above is generally believed to describe the transport mechanism of membrane transporters.30 This mechanism, however, can be accommodated through various modes of conformational changes that vary from one transporter family to another. Solcan et al.9 has suggested a “rocker-switch” mechanism2,5 for POTs, similar to those previously proposed for other MFS members.2,5,31 According to this mechanism, the entire N- and C-bundle domains of MFS transporter adopt a rocker-switch type of movement along an axis that runs along the N-,Cbundle interface.2 The rocker-switch mechanism thus can be thought of as a highly concerted global conformational change, where the N- and C-bundles of helices (each typically consisting of six helices) rotate relative to one another such that the cyto- and periplasmic (or intra- and extracellular) gating processes are negatively correlated. Recent computational studies of a MFS member, glycerol 3-phosphate transporter (GlpT), support this hypothesis and illustrate the workings of this mechanism at a detailed molecular level.31 We note that while the alternating access model has been generally accepted as the common mechanism of transporters, this mechanism does not provide specific information regarding the structural changes occurring during the alternation between the IF and OF states. Along with the rocker-switch model, other mechanisms have been proposed for various transporter families such as32 gated-pore,33,34 elevator-like,35 and toppling36 mechanisms. The gated-pore mechanism, for instance, has been proposed for certain ion-coupled transporters such as LeuT,34 which is based on the sequential use of the two gates present on the two sides of the substrate binding site.32−34 This is in contrast to the rocker-switch mechanism, which involves a concerted swivelling motion around the binding site hinge region.32 All resolved structures of GkPOT have a very similar overall conformation of the IF open state29 and thus cannot be used to elucidate the IF ↔ OF conformational transition. Doki et al.29 performed MD simulations to determine the transition mechanism of the GkPOT and investigate its conformational dynamics. On the basis of 200 ns long MD simulations of E310-unprotonated and -protonated GkPOT, they suggested that the protein is stable in its IF open conformation when the E310 residue is protonated, while the protein is in equilibrium

between the IF open and occluded forms in the absence of proton at E310 (where the E310 residue forms a salt bridge with the R43 residue).29 In addition, a 200 ns long MD simulation of E310-protonated GkPOT in complex with a phenylalanine dipeptide (FF) apparently illustrates the destabilization of the IF open state and its conversion toward an occluded conformation.29 In another study based on a series of biophysical experiments and short MD simulations, a salt bridge swapping between R36-E35 and R36-E32 was proposed to play a crucial role in structural transitions of GkPOT.37 They suggest that the R36-E32 salt bridge prefers to adopt an IF conformation, while an occluded or OF conformation is favored by the R36-E35 salt bridge.37 The reported MD simulations of POTs generally claim to be consistent with the picture presented schematically in Figure 1, where the symport mechanism is explained by the requirement of either none or both proton(s) and substrate for the accessibility of the IF ↔ OF transition. We note that, however, none of the MD simulation studies of POTs29,37,38 have captured the OF conformation and much of the arguments based on these simulations are speculative and concerned with the conformational fluctuations that are quite prevalent in MD simulations. In this study, we have investigated whether or not different simulation conditions impact the conformational dynamics of the GkPOT within the time scales accessible to typical MD simulations. To this effect, we performed all-atom MD simulations of GkPOT in an exhaustive manner. We carried out 16 different simulations of GkPOT (each for 400 ns) under different conditions in terms of the protonation state of the E310 residue and the presence/absence of substrates in the peptide binding site. The substrate binding simulations involve not only FF but also alanine dipeptide (AA), which is more similar to the substrate present in the crystal structure of GkPOT (i.e., alafosfalin),29 as well as glutamate dipeptide (EE), which is a low-affinity substrate.24 A thorough analysis of our MD simulation trajectories reveals a distinct behavior of local conformational dynamics (in terms of the substrate and binding site) with simulation conditions. Although the AA and FF substrates seem comfortable in the binding site irrespective of the protonation state of the E310, we observed a significant difference in their binding conformation. The substrates adopt a somewhat horizontal conformation similar to alafosfalin in the crystal structure in the E310-protonated state, while in the E310-unprotonated state, they move down toward the cytoplasmic side and adopt a more vertical/tilted conformation. In contrast to previous claims based on shorter simulations,29 we do not observe any statistically significant difference in the global behavior of the protein between different simulation conditions. Linear regression analysis of global quantities, however, supports the hypothesis that the GkPOT transporter functions through a concerted rocker-switch mechanism, similar to those observed in other MFS transporters.31 In the following, we will first discuss the simulation details in Theoretical Methods. Results and Discussion will provide an overview of the measures and quantities used for the analysis of our trajectories, followed by a more detailed account of the results, describing both local and global conformational fluctuations. The last section is preserved for a brief summary of our findings.

2. THEORETICAL METHODS Characterizing structural transitions of membrane transporters without compromising the detailed chemical description of 3645

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

H11 (residues 435−458), and H12 (residues 463−489). Only the Cα atoms of the residues were included in calculating the principal axes. The N- and C-bundle domains were defined as TM helices H1−H6 and H7−H12, respectively.

these systems and their environments is computationally challenging mainly due to the prohibitively long time scales involved in such processes. Here we have employed brute-force all-atom MD to study the conformational transitions of GkPOT under various conditions. First, we prepared a membrane-embedded model of the wildtype apo GkPOT in the IF state from its crystal structure (PDB code 4IKV)29 using Molecular Operating Environment (MOE)39 for preparing the protein and CHARMM-GUI40,41 for building the simulation system. All the crystal structure waters were removed, appropriate protonation states for the residues at the physiological pH (7.4) were assigned using protonate3D facility in MOE, and hydrogens and other missing atoms were added using MOE. All the glutamic acid residues were unprotonated in our simulations except for E310, which was modeled as both protonated and unprotonated. Protein was placed in the 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE) lipids, solvated in a box of TIP3P42 waters and 0.15 M NaCl using CHARMM-GUI.41 This simulation system consists of one protein, 113/107 POPE lipids in the upper/lower leaflet, 16 172 TIP3P waters, and 43/ 56 sodium/chloride ions with a total of 83 767 atoms in a box of size of ∼106 × 105 × 104 Å3 simulated with periodic boundary conditions. Given the similarity of all crystal structures, we used the same model with minor modifications to (i) protonate residue E310 and/or (ii) dock a substrate to the binding site. Three substrates were used including AA, FF, and EE, while the alafosfalin-bound GkPOT crystal structure (PDB code 4IKZ)29 was used as a template to dock the dipeptides. In total, we modeled eight different systems of GkPOT: four systems with the unprotonated-E310 and four with the protonated-E310. In each case one system is in the apo form and the other three are substrate-bound. We use CHARMM36 all-atom additive force field to describe all molecules.43,44 All MD simulations were performed with NAMD 2.10-11.45 Initially, we energy-minimized each system for 10 000 steps using conjugate gradient algorithm.46 Then, we relaxed the systems by releasing the restraints in a stepwise manner (for a total of ∼1 ns) using a procedure that is explained elsewhere.40 Further, production runs were performed starting from the relaxed models of each system for 400 ns (set 1) and repeated with a different initial velocity for another 400 ns (set 2). This initial relaxation was performed in an NVT ensemble while all production runs were performed in an NPT ensemble. Simulations were carried out using a 2 fs time step at 310 K using a Langevin integrator with a damping coefficient of γ = 0.5 ps−1. The pressure was maintained at 1 atm using the Nosé−Hoover Langevin piston method.47,48 The smoothed cutoff distance for nonbonded interactions was set to 10−12 Å, and long-range electrostatic interactions were computed with the particle mesh Ewald (PME) method.49 The trajectories were collected every 0.5 ns for the analysis, which was conducted using VMD.50 2 × 8 × 400 ns/0.5 ns = 12 800 data points were thus used for the statistical analysis presented in this work. The N-,C-bundle interdomain angle as well as the interhelical angles were calculated as the angle between the third principal axes of the respective domains/helices involved. The residues of the helices considered for these calculations were H1 (residues 23−50), H2 (residues 63−90), H3 (residues 96−113), H4 (residues 120−146), H5 (residues 154−183), H6 (residues 186−210), H7 (residues 288−322), H8 (residues 337−360), H9 (residues 369−392), H10 (residues 400−428),

3. RESULTS AND DISCUSSION 3.1. Monitoring the Conformational Dynamics. As discussed in the previous section, we have performed two independent sets of unbiased MD simulations of GkPOT (sets 1 and 2), each under various conditions including E310unprotonated (UP) and -protonated (P), in the apo and substrate-bound states, where the substrate is AA, EE, or FF. In each set (i.e., sets 1 and 2) the membrane-embedded GkPOT protein was simulated under eight different conditions including UP:apo, UP:AA, UP:FF, UP:EE, P:apo, P:AA, P:FF, and P:EE. After preparing the system for each simulation, followed by a short relaxation as described in Theoretical Methods, we performed 400 ns of unbiased MD simulations, resulting in 2 × 8 × 400 ns = 6.4 μs of MD trajectories in aggregate. The Cα root-mean-square distance (RMSD) of the protein with respect to the crystal structure (PDB code 4IKV)29 stays between 1 and 2 Å during these simulations (Figure S1 in Supporting Information), indicating the transporter is stable irrespective of the simulation condition. With two independent trajectories for each condition, one may make a more reliable comparison between various conditions. Here the duplicate simulations play the role of a control experiment to determine whether or not an observed distinct behavior of the system in a given trajectory can reliably be attributed to its distinct simulation condition. We note that although we have generated an extensive set of simulation trajectories (as compared to previous simulation studies of POT transporters29,37,38), these simulations are associated with fluctuations, which could mask the distinctive behavior of the system under different conditions. Despite this caveat, several distinctive behaviors can be observed, particularly in local dynamics of the binding site region where substrate or proton binding could result in distinct conformations within a relatively short time. On the other hand, global conformational changes of GkPOT, similar to most other transporters, are expected to be associated with much longer time scales beyond those accessible to typical unbiased MD simulations. Although one may observe fluctuations in measures associated with the global conformational features, e.g., the interdomain angles or distances, such fluctuations are not necessarily relevant to the IF ↔ OF conformational transition. Unfortunately, it is somewhat common to overinterpret such fluctuations without providing a rigorous statistical analysis, a practice that is avoided in this work. We have defined a number of measures or quantities associated with both local and global conformational features of the system under study and have monitored their behavior under different simulation conditions. Figure 2 illustrates the definition of some of these quantities including N-,C-bundle interdomain angle (Figure 2A), H1,H7 interhelical angle (Figure 2B), H5,H11 interhelical angle (Figure 2B), L4,5− L10,11 interloop distance (Figure 2C), and R43-E310 salt-bridge minimum donor−acceptor distance (Figure 2D), among which the last is a measure of local dynamics and the rest are quantities associated with global conformational features. The relative orientation of the N- and C-bundle domains of protein can be quantified using the angle between the principal axes of these domains (Figure 2A). According to the rocker 3646

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

transport mechanism of GkPOT,29 due to its relevance to proton binding (i.e., to E310) as well as its role in stabilizing certain functional states (presumably, the apo, occluded state29). In our analysis, we defined the R43-E310 distance as the minimum distance between the potential hydrogen bond donors and acceptors of R43 and E310 side chain atoms, respectively (Figure 2D). Below, we will detail our analysis of both local and global conformational changes of the simulated systems under various conditions, e.g., in terms of the quantities defined in Figure 2. 3.2. Local Conformational Changes. The E310 residue is not only a potential proton binding site; it could directly be involved in formation of the functionally relevant salt bridge R43-E310.29 To determine whether or not this salt bridge is functionally relevant, Doki et al.29 conducted functional assays with R43Q mutants, which showed the impairment of the proton-driven peptide uptake while preserving the peptidedriven counterflow activity. E310Q mutation, however, abolishes both proton-driven uptake and peptide-driven counterflow activities. Combining these observations with some MD simulations, they concluded that the R43-E310 salt bridge plays a role in the transition to the occluded, apo state (Figure 1, left side of the cycle) but not to the occluded, proton- and substrate-bound state (Figure 1, right side of the cycle). E310, on the other hand, is involved not only in the formation of the salt bridge but also in proton binding. Our extensive set of MD simulations allow for the examination of this hypothesis. As discussed below, our simulations are generally consistent with the reported functional assays associated with both proton-driven uptake and peptide-driven counterflow activities29 and may provide an explanation for the 3-fold increase in the peptide-driven counterflow activity of R43Q mutant, which was not discussed by Doki et al.29 However, our results do not provide any evidence for the R43E310 salt bridge favoring the occluded conformation, as discussed in section 3.3. 3.2.1. R43-E310 Salt Bridge. In all E310-unprotonated simulations, E310 quickly forms a salt bridge with R43, which is more or less stable in most cases. Figure 3 clearly illustrates a distinct behavior in terms of the formation of the R43-E310 salt bridge, in E310-unprotonated (Figure 3A−D) and -protonated (Figure 3E−H) systems (with the exception of the EE-bound GkPOT). The protonated E310 residue cannot form a salt bridge with the R43 residue, as expected. In our simulations, however, the EE substrate behaves differently, which is consistent with its experimentally observed low binding affinity to POTs.24 The glutamate side chains of the EE substrate could interact with the R43 residue and weaken the R43-E310 salt bridge even when the E310 residue is unprotonated (Table S1). In the other E310-unprotonated simulations, however, the R43E310 salt bridge is more or less intact. We note that since the only potential proton binding site investigated here is E310, we cannot make any conclusive statement regarding the functional role of E32 and E35, which are suggested to form salt bridges with the R36 residue in the IF and OF/occluded conformations, respectively.37 In our simulations of the IF state, however, the R36 residue consistently forms a salt bridge with E35 rather than E32, while K136 forms salt bridges with both E32 and E35. This interaction pattern was consistently observed in all of the eight systems simulated as quantified using the minimum donor−acceptor distance between the saltbridge forming residues (Table S2).

Figure 2. GkPOT transporter in its IF conformation (a representative snapshot from UP:apo, set 1). The definitions of some of the measures quantifying important global (A−C) and local (D) conformational features are illustrated. (A) N- and C-bundles are colored green and orange, respectively. (B, D) TM helices are colored yellow (H1,H7), green (H2,H8), orange (H5,H11), and gray (H3,H4,H6,H9,H10,H12). (C) The L4,5 and L10,11 loops are colored blue.

switch mechanism, the relative rotational changes of the N- and C-bundles are the main mode of motion that provides the structural basis of the alternating access mechanism; thereby the N-,C-bundle interdomain angle is a useful quantity to monitor the global conformational changes associated with the IF → OF transition. In addition, interhelical angles can be monitored using the angle between the roll (i.e., the third principal) axes of the helices (Figure 2B). The periplasmic gate of GkPOT is formed by the helices H1,H2 (from N-bundle) and H7,H8 (from C-bundle). Therefore, the interhelical angles such as those between H1 and H7 (Figure 2B), or H2 and H8 are expected to reflect the global conformational changes associated with the opening/closing of the periplasmic gate. Similarly, the cytoplasmic gate is formed by the TM helices H4,H5 (from N-bundle) and H10,H11 (from C-bundle), and interhelical angles between H4 and H10, or H5 and H11 (Figure 2B) reflect the global conformational changes associated with the opening/closing of the cytoplasmic gate. These measures are similar to those recently used to describe the global conformational changes of GlpT31 and other membrane transporters.51,52 We note that along with these less conventional measures for monitoring the IF → OF transition associated conformational changes, we have monitored more conventional quantities such as RMSD (with respect to our initial model based on the crystal structure29) and root-mean-square fluctuations (RMSF) as well as certain interloop distances. For instance, the loop between H4 and H5 (L4,5) and the loop between H10 and H11 (L10,11) may approach one another during the closure of cytoplasmic gate and the L4,5−L10,11 interloop distance (Figure 2C) could be used to monitor the closure of this gate.29 In terms of the local dynamics of the binding site, the R43E310 salt-bridge is proposed to play an important role in the 3647

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

Figure 3. R43-E310 salt bridge fluctuations in GkPOT equilibrium simulations monitored by their minimum donor−acceptor distance. The left (A− D) and right panels (E−H) represent the unprotonated- and protonated-E310 systems, respectively.

3.2.3. Peptide-Driven Counterflow Increase in the R43Q Mutant. An observation based on our analysis of R43-E310 distance is that the R43-E310 salt bridge could form not only in the absence but also in the presence of the substrate. We believe this is consistent with the experimental observation that the R43Q mutation not only preserves but also increases the peptide-driven counterflow activity by a factor of ∼3.29 We hypothesize that when the proton concentration is low and a substrate is bound to the protein prior to the proton, the formation of the R43-E310 salt bridge could lower the proton affinity to the E310 residue. Therefore, the R43Q mutant would have a higher affinity to proton. In a peptide-driven process, where the substrate is more likely to bind to the protein as compared to the proton, the absence of the R43-E310 salt bridge (due to the R43Q mutation) could boost the E310 proton binding (by eliminating the competing R43 basic side chain). Our simulations are consistent with this hypothesis, since they provide evidence for the formation of R43-E310 salt bridge in the presence of the substrate. 3.2.4. Substrate Binding Conformation. While we used R43-E310 salt bridge distance as a representative quantity to describe the local conformational changes of binding site, several substrate-related quantities were measured as well, as summarized in Figure 4. For the sake of brevity, the data are

3.2.2. EE Unbinding Event. We note that while in the UP:EE set 2 simulation (Figure 3D) the R43-E310 salt bridge is quickly disrupted by the presence of the EE substrate (due to the interaction of EE and R43 side chains), the EE substrate is not able to disrupt this salt bridge in the UP:EE set 1 simulation. This is due to the fact that the EE substrate is not stably bound to the protein in the UP:EE set 1 simulation and leaves the binding site within the first 50 ns of these unbiased simulations. The EE peptide is experimentally known to have a low affinity to GkPOT homolog PepTSt(IC50 ≈ 0.6 mM) as compared to the FF (0.005 mM) and AA (0.04 mM) peptides.24 Therefore, it is not surprising that the only observed unbinding in these 16 simulations was associated with the EE. The unbinding of the EE substrate, however, does not completely rescue the R43-E310 salt bridge from breaking. Although this system should eventually behave similarly to the UP:apo system (Figure 3A), the presence and subsequently the unbinding of the substrate have already hydrated the binding site to an extent that the R43-E310 salt bridge can be disrupted for ∼100 ns before it forms again. This observation indicates that even local conformational changes could be slow and may not be fully understood based on short simulations, typically used in unbiased MD studies. 3648

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

Figure 4. Substrate binding conformation described in terms of the average (squares) and standard deviation (error bars) of various measures.

presented in terms of the average and standard deviation (shown as error bars) associated with each quantity as obtained from individual trajectories. Figure 4A and Figure 4B compare the substrate heavy-atom RMSD obtained from various simulations. The overall RMSD is obtained by aligning the protein with respect to the crystal structure and calculating the substrate RMSD with respect to its initial conformation in the binding site (modeled using the alafosfalin in the crystal structure (PDB code 4IKZ)29). The internal RMSD is calculated by aligning the substrate with its own initial conformation and calculating its RMSD with respect to this conformation. Therefore, the substrate overall RMSD represents a combination of its internal conformational changes with respect to the initial model and its overall rotational and translational movement within the binding site. The AA overall conformation seems to be more stable in the E310-protonated simulations than in the -unprotonated ones, which may indicate the substrate is more comfortable in the presence of the proton in the binding site. However, this trend is not consistently

observed in all FF and EE simulations (e.g., UP:FF, set 1 as well as both EE, Set 2 simulations). While the RMSD based analysis of the substrate dynamics is not conclusive, we can observe more distinctive behavior for different conditions when the analysis is more specific. For instance, the Z position of the substrate mass center shows an interestingly distinct behavior in E310-unprotonated and -protonated simulations (Figure 4C). The substrate is placed closer to the periplasm in the binding site by about 2−4 Å irrespective of the type of substrate used (excluding the EE substrate that leaves the binding site). In terms of the internal substrate conformation, monitoring the torsion angles reveal that the substrate stays within the extended/β-sheet/PPII region of the Ramachandran map with the exception of the EE peptide that stays bound to the E310-unprotonated GkPOT (see Ψ1 torsion angle values shown in Figure 4D). This may indicate that unlike the AA and FF dipeptides, the EE is not comfortable in the binding site in its extended conformation, when E310 is unprotonated, further justifying the unbinding of the extended EE in the UP:EE, set 1 simulation. Note that the 3649

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

Figure 5. Substrate conformation in UP:AA (A) and P:AA (B) systems, illustrated by representative snapshots. Residues with highest interaction probability with the substrate are highlighted, and the R43-E310 salt bridge is shown. The salt bridge residues in panel A are faded, since they are not interacting with the substrate. The substrate (AA) is represented by the ball-and-stick model. TM helices H1,H7 and H5,H11 are colored yellow and orange, respectively.

of the binding site and the protein based on short MD simulations, as discussed in section 3.3. Table S1 provides information on the substrate binding site and the residues interacting with the substrate. While the prominent role of E413 residue (with the exception of EE simulations) is clear, several other residues are involved such as N342 and Y40. When E310 is protonated, both R43 and E310 interact with the substrate. We note that while the substrate binding conformations shown in Figure 5 are representative of the conformations observed in our simulations, we cannot exclude the possibility of other binding conformations, which could potentially be more stable on longer time scales. In our simulations, we did observe other short-lived binding conformations, some of which are shown in Figure S3; however, the time scales associated with the function of GkPOT is much longer than those observed in our MD studies. 3.3. Global Conformational Changes. In order to quantify the global conformational changes of protein, particularly those associated with the IF → OF transition, we have monitored the interdomain orientational changes of a number of helical bundles or individual helices such as those defined in Figure 2A−C. We note that the H1,H7 and H5,H11 interhelical angles here represent the intuitive measures associated with peri- and cytoplasmic gating. However, other interhelical angles involving TM helices H3, H6, H9, and H12 could indirectly be relevant to the rocker-switch mechanism. The conventional measures such as RMSD (Figure S1) were also used in our analysis. However, the more specific measures such as those described above provide a more intuitive description of the global conformational changes that are relevant to the state transition. We have also included the L4,5− L10,11 interloop distance (Figure 2C) in our analysis as a measure for cytoplasmic gating as used by Doki et al.,29 primarily to compare the results of the two MD studies. For a more direct measure of gating behavior, we have also quantified the water accessibility of both cyto- and periplasmic gates (Figure 6E,F), which could be a measure of substrate accessibility. This analysis is particularly important when the IF → OF conformational transition is concerned. The exact definition of the peri-/cytoplasmic water count is illustrated in Figure S4, where the water content of the individual simulations is shown as histograms along the pore. 3.3.1. Comparing Global Conformational Changes. Our results do not show any statistically significant difference in the global behavior of protein under different conditions (Figure

internal RMSD analysis does not reveal the distinct conformation of the substrate in the UP:EE, set 1 simulation, due to the typical degeneracies associated with the singlereference RMSDs. 3.2.5. Substrate Binding Orientation. The orientation of the substrate in the binding site was quantified using the tilt and azimuth angles of the substrate with respect to a fixed reference frame aligned with the protein. The tilt angle quantifies the angle between the membrane normal and the vector connecting the two Cα atoms of the substrate, with positive/negative values indicating that the C-terminal of the peptide is directed toward the peri-/cytoplasm. In most cases, the substrate is tilted toward the cytoplasm (as defined above) as shown in Figure 4E. The exception is the substrate in P:AA and P:FF simulations, where it fluctuates around a lateral orientation. This is another indication of the existence of two substantially different binding conformations for E310-protonated and -unprotonated systems (with the usual exception of the EE peptide). The picture that emerges from these analyses is that the substrate (excluding the EE) is placed toward the lower part of the binding site (i.e., toward the cytoplasm) and tilted downward in a vertical conformation when the E310 is unprotonated (Figure 5A). Protonation of the E310, however, makes the substrate uncomfortable in this conformation, while breaking of the R43-E310 salt bridge allows for its upward movement. The resulting conformation, however, is almost lateral (Figure 5B) with considerable fluctuations in the tilt angle. These fluctuations are not necessarily indicative of weak binding, as the azimuth angle does not reflect any instability (Figure 4F). 3.2.6. Substrate Binding Site. If we consider the number of residues that are interacting with the substrate with high frequencies (e.g., greater than 70% as shown in Table S1), one may see a clear trend, where the E310-protonated systems are in stable contact with more residues as compared to the E310unprotonated systems. This may indicate a more favorable interaction with the substrate when the E310 residue is protonated. We note that the accurate quantification of the binding affinities requires employing enhanced sampling methods. A common hypothesis with regard to symporters such as GkPOT is that the co-transported species (here the proton) lock the protein conformation in an (IF or OF) open state, thereby promoting the substrate binding (Figure 1, middle). It is difficult, however, to speculate on the interdependence of the local and global conformational changes 3650

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

Figure 6. GkPOT global conformational behavior described in terms of the average (squares) and standard deviation (error bars) of various measures including N-,C-bundle interdomain angle (A), H1,H7 interhelical angle (B), H5,H11 interhelical angle (C), L4,5−L10,11 distance (D), and water count at the cytoplasmic (E) and periplasmic (F) gates. Green arrows represent the values associated with the crystal structure (PDB code 4IKV).

6). This is indeed expected, since the time scale of the IF → OF transition is much longer than those currently accessible to typical unbiased simulations such as those presented here. Figure 6 clearly shows that neither protonation of the E310 residue nor the binding of a substrate significantly changes the protein global conformational dynamics within the time scales of our simulations. 3.3.2. Comparing Individual Trajectories. We note that one may see differences between individual trajectories if “select” trajectories are considered. However, only the repeated simulations can provide enough statistics to make a reliable comparison between different conditions. Figure 7A and Figure 7B show examples of the apparently different behaviors of GkPOT under two different conditions, when only set 1 is

taken into consideration. Set 2 simulations, however, do not show any meaningful difference in the behavior of these quantities (Figure 7C and Figure 7D). Similarly, one can see a similar behavior for other global quantities such as interhelical angles (Figure S5). Obviously, we do not claim that GkPOT behaves the same way under various conditions. We simply point out the fact that one cannot make a connection between the stochastic, short-time scale behavior of the system and its functionally relevant, long-time scale behavior unless the difference in fluctuations is statistically significant. We note that Doki et al.29 have made observations with regard to the different global behavior of the GkPOT protein under various conditions (e.g., L4,5−L10,11 interloop distance). These observations are not inconsistent with ours; however, the observed 3651

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

Figure 7. GkPOT global conformational dynamics described by an interdomain angle (A, C) and an interloop distance (B, D) for select simulations.

3.3.4. Indication of a Concerted Motion. The TM helices H6 and H12 are placed quite apart from each other and are not directly involved in the cyto- or periplasmic gating. Their interhelical angle, however, fluctuates significantly similar to that of H1,H7 or H5,H11 discussed above and is correlated with both quantities with a correlation coefficient of r ≈ 0.7 and −0.7, respectively (Figure 8B,C). This observation indicates that the concerted interhelical fluctuations are not limited to the helices that are directly involved in the gating. The picture that emerges is consistent with a global, concerted motion in line with the rocker-switch mechanism. The N-,C-bundle interdomain angle also strongly correlates with the H1,H7 interhelical angle with a correlation coefficient of r ≈ 0.85 (Figure 8D). We note that this particular correlation is partly inherent, due to the fact that H1 and H7 helices are part of the N- and C-bundles, respectively. However, the observed correlations between interhelical angles (Figure 8A−C) imply that the strong correlation seen in Figure 8D is not solely due to this inherent correlation. One may think of the deviation of the absolute value of these correlation coefficients from 1.0 as a measure of the deviation from an idealized rocker-switch mechanism, where all helices participate in a “concerted motion” to result in a global transition between the IF and OF states. The conformational changes that occur in our simulations, however, are only partly related to this concerted motion. For instance, H1,H7 interhelical angle fluctuations explain only about 30% of the total fluctuations as quantified by coefficient of determination, r2, where r is the correlation coefficient between the H1,H7 interhelical angle and Cα RMSD of the protein (Figure 8E). This particular correlation is also inherent, since the RMSD is generally correlated with any other measure. The correlation between RMSD and different measures is particularly interesting, since it can be used to quantify the contribution of fluctuations associated with those

differences are not statistically significant based on the more extensive simulations reported here (2 × 400 ns per system). Doki et al. make other observations regarding the distinctive behavior of GkPOT under different conditions in terms of RMSF of the L4,5 region (or the H4−H5 hairpin) and the angle between the H4,H5 helices and the C-bundle.29 As illustrated in Figure S6, we do not observe any significantly distinct behavior for different simulation conditions in terms of these measures, considering the extensive set of simulations conducted here. 3.3.3. Indication of the Rocker-Switch Mechanism. Although no significant global conformational change is observed in our simulations that can be directly indicative of a state transition, significant correlations are observed between certain global quantities, which could be relevant to the alternating access and rocker-switch mechanisms (Figure 8). For instance, one may observe a statistically significant correlation between quantities such as H1,H7 interhelical angle (which is directly involved in the periplasmic gating) and H5,H11 interhelical angle (which is directly involved in the cytoplasmic gating) as illustrated in Figure 8A. A linear correlation of r ≈ −0.76 (where r is the linear correlation coefficient) between the two quantities indicates that the cytoand periplasmic gating processes are anticorrelated, which is consistent with both alternating access and rocker switch mechanisms. The correlations between several quantities associated with the global conformational changes are illustrated in Figure 8, where the estimated correlation coefficients are shown as well. Table S3 summarizes the correlation coefficients of a select number of quantities representing global conformational fluctuations, while the correlation coefficients obtained from individual simulation trajectories are shown in Figure S7. 3652

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

Figure 8. Correlations between various measures associated with global conformational changes of GkPOT. Each simulation system is represented by a color (see the color bar on the right).

3.3.5. Correlation between Local and Global Fluctuations. While several measures associated with the global conformational fluctuations correlate strongly with each other, certain quantities associated with the local fluctuations of the substrate and the binding site show some degree of correlation as well. Table S4 shows examples of the correlation coefficients between different local quantities. For instance, the Z position of the substrate weakly correlates with its tilt angle with a correlation coefficient of 0.46. However, the stronger local correlations are often trivial. For instance, both the Z position of the substrate and its tilt angle correlate strongly with its overall RMSD, as expected. We did not find any statistically meaningful correlation between any local and global quantities measured. The largest r2 value between any local and global quantity was about 0.04 (r ≈ 0.2), which is almost negligible.

measures to the variance using the coefficient of determination. The relatively weak correlation observed between the RMSD and the H1,H7 interhelical angle is partly due to the inconsistent behavior of different trajectories. For instance, the two sets of UP:apo simulations are associated with correlation coefficients 0.44 and −0.76 (for sets 1 and 2, respectively) between the two quantities. The other 14 simulations have values ranging from approximately −0.08 to ∼0.78 for this correlation (with r2 values ranging from ∼1% to ∼61%). Limited statistics could thus be quite misleading regarding the contributions of various TM helical fluctuations to the variance. Finally, Figure 8F reveals the strong correlation between the L4,5−L10,11 interloop distance and the H1,H7 interhelical angle as mentioned above. 3653

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

respectively. The substrate binding conformation, however, is quite different in the two conditions. The substrate binds to the E310-protonated GkPOT in a lateral manner similar to the alafosfalin captured in the crystal structure.29 When the E310 residue is unprotonated, however, the substrate is moved toward the cytoplasmic side and is tilted with respect to the membrane normal with its C terminal pointing toward the cytoplasm. This behavior is consistently observed not only in both independent sets of simulations but also for both AA and FF substrates. Unlike the local conformational dynamics of GkPOT, we do not observe any statistically significant distinctive behavior in terms of its global conformational changes. However, linear regression analysis of quantities associated with the global fluctuations gives rise to a general picture of a concerted motion consistent with the rocker-switch mechanism of other MFS transporters.2,5,31 Combining all 16 independent sets of simulations (i.e., 6.4 μs), we can observe a statistically meaningful anticorrelation between the fluctuations of the key transmembrane helices involved in the cytoplasmic gating (quantified by the H1,H7 interhelical angle) and those involved in the periplasmic gating (quantified by the H5,H11 interhelical angle). These correlations are not limited to the helices that are directly involved in the gating. These correlations support the possibility of a concerted mode of motion that involves other transmembrane helices such as H6 and H12. Similarly, meaningful but weaker correlations are observed between certain local variables (e.g., the Z position of substrate and its tilt angle). Interestingly, however, we do not observe any significant correlations between local and global variables (with largest correlation coefficient being ∼0.2). The presented study provides a detailed picture of the substrate binding conformation and the local effect of the proton binding, both in terms of the salt bridge networks and the substrate binding. However, the time scales of these simulations are much shorter than those associated with the IF ↔ OF conformational transition. These simulations could provide the basis for the application of enhanced sampling techniques, similar to those recently conducted for the study of global transitions in another MFS transporter.31

These observations call into question the implicit underlying assumption of many MD studies that the global conformational changes of proteins are triggered by local binding events on nanosecond time scales. The autocorrelation functions associated with both local and global quantities as shown in Figure S8 indicate that these quantities are associated with autocorrelation times of ∼50 to 300 ns. Therefore, much longer simulations are required for accurate statistical analysis of both local and global conformational changes in these systems. Recent MD simulations of various membrane transport proteins have attempted to provide a uniquely detailed dynamic picture of membrane transport at an atomic level53−63 because of the ever increasing power of computers, improved force fields,43,44,64−68 and high-throughput modeling tools such as CHARMM-GUI.69−73 However, the accurate description of global conformational changes based on short MD simulations has proved challenging. What we have reported in this work sheds light on the distinct nature of characterizing local and global conformational dynamics using typical MD simulations of membrane transport proteins.

4. CONCLUSION A fundamental aspect of the transport mechanism in membrane transporters is the coupling between local conformational changes induced by substrate or ion binding and protein global conformational changes such as IF ↔ OF transition. Many MD studies are aimed at describing these couplings using simulations of the protein under different conditions, e.g., in the presence or absence of a substrate or ion. While the current study follows a similar approach, it presents a critical view of deciphering information regarding the global conformational changes based on local conformational dynamics. Using the most extensive set of all-atom MD simulations of a POT transporter to date, we have demonstrated that the distinct local dynamics of the binding site and the substrate under various conditions does not result in a statistically meaningful distinctive behavior at the global level of protein conformational dynamics. The global conformational variations, on the other hand, are primarily due to the thermal fluctuations present across the simulations irrespective of the simulation conditions. The 16 independent MD trajectories of membraneembedded GkPOT generated under eight different conditions allowed us to study the effect of substrate/proton binding on both global and local conformational changes and determine whether or not such simulations could be used to make conclusive statements with respect to either type of conformational changes. In these simulations, the proton binding site residue E310 was simulated in both its protonated and unprotonated states, either with or without a substrate. Three different substrates were used in the substrate-bound simulations including the AA (which is quite similar to alafosfalin captured in the GkPOT crystal structure29) as well as FF and EE, representing a high- and a low-affinity dipeptide, respectively. A careful analysis of the local conformational dynamics in terms of both the substrate and the binding site reveals a clearly distinct behavior in the absence and presence of the proton at the E310 side chain. Unlike the EE that represents a low-affinity substrate, the other substrates seem to be comfortable in both E310-protonated and -unprotonated conditions. This is despite the fact that the R43-E310 salt bridge is consistently broken/ formed in the E310-protonated/-unprotonated condition,



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b09733. Figures S1−S8 and Tables S1−S4 providing additional analysis based on our MD trajectories (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +1 (479) 575-6459. ORCID

Mahmoud Moradi: 0000-0002-0601-402X Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This research is supported by the University of Arkansas, Fayetteville. This research is part of the Blue Waters sustainedpetascale computing project, which is supported by the National Science Foundation (Awards OCI-0725070 and 3654

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

(20) Lyons, J. A.; Parker, J. L.; Solcan, N.; Brinth, A.; Li, D.; Shah, S. T.; Caffrey, M.; Newstead, S. Structural basis for polyspecificity in the POT family of proton-coupled oligopeptide transporters. EMBO Rep. 2014, 15, 886−893. (21) Jensen, J. M.; Simonsen, F. C.; Mastali, A.; Hald, H.; Lillebro, I.; Diness, F.; Olsen, L.; Mirza, O. Biophysical characterization of the proton-coupled oligopeptide transporter YjdL. Peptides 2012, 38, 89− 93. (22) Prabhala, B. K.; Aduri, N. G.; Hald, H.; Mirza, O. Investigation of the substrate specificity of the proton coupled peptide transporter PepTSo from Shewanella oneidensis. Int. J. Pept. Res. Ther. 2015, 21, 1− 6. (23) Guettou, F.; Quistgaard, E. M.; Tresaugues, L.; Moberg, P.; Jegerschold, C.; Zhu, L.; Jong, A. J.; Nordlund, P.; Loew, C. Structural insights into substrate recognition in proton-dependent oligopeptide transporters. EMBO Rep. 2013, 14, 804−810. (24) Newstead, S. Molecular insights into proton coupled peptide transport in the PTR family of oligopeptide transporters. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850, 488−499. (25) Smith, D. E.; Clémençon, B.; Hediger, M. A. Proton-coupled oligopeptide transporter family SLC15: Physiological, pharmacological and pathological implications. Mol. Aspects Med. 2013, 34, 323−336. (26) Meredith, D.; Boyd, C. A. R. Structure and function of eukaryotic peptide transporters. Cell. Mol. Life Sci. 2000, 57, 754−778. (27) Rautio, J.; Kumpulainen, H.; Heimbach, T.; Oliyai, R.; Oh, D.; Jarvinen, T.; Savolainen, J. Prodrugs: design and clinical applications. Nat. Rev. Drug Discovery 2008, 7, 255−270. (28) Brandsch, M. Drug transport via the intestinal peptide transporter PepT1. Curr. Opin. Pharmacol. 2013, 13, 881−887. (29) Doki, S.; Kato, H. E.; Solcan, N.; Iwaki, M.; Koyama, M.; Hattori, M.; Iwase, N.; Tsukazaki, T.; Sugita, Y.; Kandori, H.; et al. Structural basis for dynamic mechanism of proton-coupled symport by the peptide transporter POT. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 11343−11348. (30) Jardetzky, O. Simple allosteric model for membrane pumps. Nature 1966, 211, 969−970. (31) Moradi, M.; Enkavi, G.; Tajkhorshid, E. Atomic-level characterization of transport cycle thermodynamics in the glycerol-3phosphate:phosphate antiporter. Nat. Commun. 2015, 6, 8393. (32) Slotboom, D. J. Structural and mechanistic insights into prokaryotic energy-coupling factor transporters. Nat. Rev. Microbiol. 2014, 12, 79−87. (33) Korkhov, V. M.; Mireku, S. A.; Locher, K. P. Structure of AMPPNP-bound vitamin B12 transporter BtuCD-F. Nature 2012, 490, 367−372. (34) Forrest, L. R.; Rudnick, G. The rocking bundle: A mechanism for ion-coupled solute flux by symmetrical transporters. Physiology 2009, 24, 377−386. (35) Ryan, R.; Vandenberg, R. Elevating the alternating-access model. Nat. Struct. Mol. Biol. 2016, 23, 187−189. (36) Xu, K.; Zhang, M.; Zhao, Q.; Yu, F.; Guo, H.; Wang, C.; He, F.; Ding, J.; Zhang, P. Crystal structure of a folate energy-coupling factor transporter from Lactobacillus brevis. Nature 2013, 497, 268−271. (37) Aduri, N. G.; Prabhala, B. K.; Ernst, H. A.; Jørgensen, F. S.; Olsen, L.; Mirza, O. Salt bridge swapping in the EXXERFXYY motif of proton-coupled oligopeptide transporters. J. Biol. Chem. 2015, 290, 29931−29940. (38) Fowler, P. W.; Orwick-Rydmark, M.; Radestock, S.; Solcan, N.; Dijkman, P. M.; Lyons, J. A.; Kwok, J.; Caffrey, M.; Watts, A.; Forrest, L.; et al. Gating topology of the proton-coupled oligopeptide symporters. Structure 2015, 23, 290−301. (39) Molecular Operating Environment (MOE), version 2013.08; Chemical Computing Group Inc.: Montreal, QC, Canada, 2016. (40) Jo, S.; Kim, T.; Im, W. Automated builder and database of protein/membrane complexes for molecular dynamics simulations. PLoS One 2007, 2, e880. (41) Lee, J.; Cheng, X.; Swails, J. M.; Yeom, M. S.; Eastman, P. K.; Lemkul, J. A.; Wei, S.; Buckner, J.; Jeong, J. C.; Qi, Y.; et al. CHARMM-GUI input generator for NAMD, GROMACS, AMBER,

ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at UrbanaChampaign and its National Center for Supercomputing Applications. This research is also supported by the Arkansas High Performance Computing Center which is funded through multiple National Science Foundation grants and the Arkansas Economic Development Commission.



REFERENCES

(1) Marger, M.; Saier, M. H., Jr. A major superfamily of transmembrane facilitators that catalyse uniport, symport and antiport. Trends Biochem. Sci. 1993, 18, 13−20. (2) Law, C. J.; Maloney, P. C.; Wang, D.-N. Ins and outs of major facilitator superfamily antiporters. Annu. Rev. Microbiol. 2008, 62, 289− 305. (3) Yan, N. Structural advances for the major facilitator superfamily (MFS) transporters. Trends Biochem. Sci. 2013, 38, 151−159. (4) Abramson, J.; Smirnova, I.; Kasho, V.; Verner, G.; Kaback, H. R.; Iwata, S. Structure and mechanism of the lactose permease of Escherichia coli. Science 2003, 301, 610−615. (5) Huang, Y.; Lemieux, M. J.; Song, J.; Auer, M.; Wang, D.-N. Structure and mechanism of the glycerol-3-phosphate transporter from Escherichia coli. Science 2003, 301, 616−620. (6) Yin, Y.; He, X.; Szewczyk, P.; Nguyen, T.; Chang, G. Structure of the multidrug transporter EmrD from Escherichia coli. Science 2006, 312, 741−744. (7) Dang, S.; Sun, L.; Huang, Y.; Lu, F.; Liu, Y.; Gong, H.; Wang, J.; Yan, N. Structure of a fucose transporter in an outward-open conformation. Nature 2010, 467, 734−738. (8) Newstead, S.; Drew, D.; Cameron, A. D.; Postis, V. L.; Xia, X.; Fowler, P. W.; Ingram, J. C.; Carpenter, E. P.; Sansom, M. S.; McPherson, M. J.; et al. Crystal structure of a prokaryotic homologue of the mammalian oligopeptide-proton symporters, PepT1 and PepT2. EMBO J. 2011, 30, 417−426. (9) Solcan, N.; Kwok, J.; Fowler, P. W.; Cameron, A. D.; Drew, D.; Iwata, S.; Newstead, S. Alternating access mechanism in the POT family of oligopeptide transporters. EMBO J. 2012, 31, 3411−3421. (10) Sun, L.; Zeng, X.; Yan, C.; Sun, X.; Gong, X.; Rao, Y.; Yan, N. Crystal structure of a bacterial homologue of glucose transporters GLUT1−4. Nature 2012, 490, 361−366. (11) Pedersen, B. P.; Kumar, H.; Waight, A. B.; Risenmay, A. J.; RoeZurz, Z.; Chau, B. H.; Schlessinger, A.; Bonomi, M.; Harries, W.; Sali, A.; et al. Crystal structure of a eukaryotic phosphate transporter. Nature 2013, 496, 533−536. (12) Quistgaard, E. M.; Löw, C.; Moberg, P.; Trésaugues, L.; Nordlund, P. Structural basis for substrate transport in GLUThomology family of monosaccharide transporters. Nat. Struct. Mol. Biol. 2013, 20, 766−768. (13) Iancu, C. V.; Zamoon, J.; Woo, S. B.; Aleshin, A.; Choe, J.-y. Crystal structure of a glucose/H+ symporter and its mechanism of action. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 17862−17867. (14) Parker, J. L.; Newstead, S. Molecular basis of nitrate uptake by the plant nitrate transporter NRT1.1. Nature 2014, 507, 68−72. (15) Sun, J.; Bankston, J. R.; Payandeh, J.; Hinds, T. R.; Zagotta, W. N.; Zheng, N. Crystal structure of the plant dual-affinity nitrate transporter NRT1.1. Nature 2014, 507, 73−77. (16) Kumar, H.; Kasho, V.; Smirnova, I.; Finer-Moore, J. S.; Kaback, H. R.; Stroud, R. M. Structure of sugar-bound LacY. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 1784−1788. (17) Wisedchaisri, G.; Park, M. S.; Iadanza, M. G.; Zheng, H.; Gonen, T. Proton-coupled sugar transport in the prototypical major facilitator superfamily protein XylE. Nat. Commun. 2014, 5, 4521. (18) Paulsen, I. T.; Skurray, R. A. The POT family of transport proteins. Trends Biochem. Sci. 1994, 19, 404. (19) Ito, K.; Hikida, A.; Kawai, S.; Lan, V. T. T.; Motoyama, T.; Kitagawa, S.; Yoshikawa, Y.; Kato, R.; Kawarasaki, Y. Analysing the substrate multispecificity of a proton-coupled oligopeptide transporter using a dipeptide library. Nat. Commun. 2013, 4, 2502. 3655

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656

Article

The Journal of Physical Chemistry B

(61) Mondal, S.; Khelashvili, G.; Shi, L.; Weinstein, H. The cost of living in the membrane: a case study of hydrophobic mismatch for the multi-segment protein LeuT. Chem. Phys. Lipids 2013, 169, 27−38. (62) Wen, P.-C.; Verhalen, B.; Wilkens, S.; Mchaourab, H. S.; Tajkhorshid, E. On the origin of large flexibility of P-glycoprotein in the inward-facing state. J. Biol. Chem. 2013, 288, 19211−19220. (63) Cheng, M. H.; Bahar, I. Complete mapping of substrate translocation highlights the role of LeuT N-terminal segment in regulating transport cycle. PLoS Comput. Biol. 2014, 10, e1003879. (64) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671−690. (65) Yesselman, J. D.; Price, D. J.; Knight, J. L.; Brooks, C. L., III MATCH: An atom−typing toolset for molecular mechanics force fields. J. Comput. Chem. 2012, 33, 189−202. (66) Yu, W.; He, X.; Vanommeslaeghe, K.; MacKerell, A. D., Jr. Extension of the CHARMM general force field to sulfonyl-containing compounds and its utility in biomolecular simulations. J. Comput. Chem. 2012, 33, 2451−2468. (67) Burger, S. K.; Ayers, P. W.; Schofield, J. Efficient parameterization of torsional terms for force fields. J. Comput. Chem. 2014, 35, 1438−1445. (68) Coimbra, J. T.; Sousa, S. F.; Fernandes, P. A.; Rangel, M.; Ramos, M. J. Biomembrane simulations of 12 lipid types using the general amber force field in a tensionless ensemble. J. Biomol. Struct. Dyn. 2014, 32, 88−103. (69) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A webbased graphical user interface for CHARMM. J. Comput. Chem. 2008, 29, 1859−1865. (70) Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W. CHARMM-GUI membrane builder for mixed bilayers and its application to yeast membranes. Biophys. J. 2009, 97, 50−58. (71) Wu, E.; Cheng, X.; Jo, S.; Rui, H.; Song, K.; Dávila-Contreras, E.; Qi, Y.; Lee, J.; Monje-Galvan, V.; Venable, R.; et al. CHARMMGUI membrane builder toward realistic biological membrane simulations. J. Comput. Chem. 2014, 35, 1997−2004. (72) Qi, Y.; Cheng, X.; Han, W.; Jo, S.; Schulten, K.; Im, W. CHARMM-GUI PACE CG builder for solution, micelle, and bilayer coarse-grained simulations. J. Chem. Inf. Model. 2014, 54, 1003−1009. (73) Qi, Y.; Ingólfsson, H. I.; Cheng, X.; Lee, J.; Marrink, S. J.; Im, W. CHARMM-GUI Martini maker for coarse-grained simulations with the Martini force field. J. Chem. Theory Comput. 2015, 11, 4486−4494.

OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. J. Chem. Theory Comput. 2016, 12, 405−413. (42) 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. (43) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr.; Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: Validation on six lipid types. J. Phys. Chem. B 2010, 114, 7830−7843. (44) 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 dngles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (45) 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. (46) Reid, J. K. In Large Sparse Sets of Linear Equations; Reid, J. K., Ed.; Academic Press: London, 1971; pp 231−254. (47) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177−4189. (48) 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. (49) Darden, T.; York, D.; Pedersen, L. G. Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (50) Humphrey, W.; Dalke, A.; Schulten, K. VMD − visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (51) Moradi, M.; Tajkhorshid, E. Mechanistic picture for conformational transition of a membrane transporter at atomic resolution. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 18916−18921. (52) Moradi, M.; Tajkhorshid, E. Computational recipe for efficient description of large-scale conformational changes in biomolecular systems. J. Chem. Theory Comput. 2014, 10, 2866−2880. (53) Khafizov, K.; Perez, C.; Koshy, C.; Quick, M.; Fendler, K.; Ziegler, C.; Forrest, L. R. Investigation of the sodium-binding sites in the sodium-coupled betaine transporter BetP. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, E3035−E3044. (54) Vanni, S.; Campomanes, P.; Marcia, M.; Rothlisberger, U. Ion binding and internal hydration in the multidrug resistance secondary active transporter NorM investigated by molecular dynamics simulations. Biochemistry 2012, 51, 1281−1287. (55) Mehmood, S.; Domene, C.; Forest, E.; Jault, J.-M. Dynamics of a bacterial multidrug ABC transporter in the inward- and outwardfacing conformations. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 10832− 10836. (56) 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. (57) Zomot, E.; Bahar, I. A Conformational switch in a partially unwound helix selectively determines the pathway for substrate release from the carnitine/γ-butyrobetaine antiporter CaiT. J. Biol. Chem. 2012, 287, 31823−31832. (58) Bisha, I.; Laio, A.; Magistrato, A.; Giorgetti, A.; Sgrignani, J. A candidate ion-retaining state in the inward-facing conformation of sodium/galactose symporter: Clues from atomistic simulations. J. Chem. Theory Comput. 2013, 9, 1240−1246. (59) Zhao, C.; Noskov, S. Y. The molecular mechanism of iondependent gating in secondary transporters. PLoS Comput. Biol. 2013, 9, e1003296. (60) Heinzelmann, G.; Bastug, T.; Kuyucak, S. Mechanism and energetics of ligand release in the aspartate transporter GltPh. J. Phys. Chem. B 2013, 117, 5486−5496. 3656

DOI: 10.1021/acs.jpcb.6b09733 J. Phys. Chem. B 2017, 121, 3644−3656