Visualizing functional motions of membrane transporters with molecular dynamics simulations.

Computational modeling and molecular simulation techniques have become an integral part of modern molecular research. Various areas of molecular sciences continue to benefit from, indeed rely on, the unparalleled spatial and temporal resolutions offered by these technologies, to provide a more complete picture of the molecular problems at hand. Because of the continuous development of more efficient algorithms harvesting ever-expanding computational resources, and the emergence of more advanced and novel theories and methodologies, the scope of computational studies has expanded significantly over the past decade, now including much larger molecular systems and far more complex molecular phenomena. Among the various computer modeling techniques, the application of molecular dynamics (MD) simulation and related techniques has particularly drawn attention in biomolecular research, because of the ability of the method to describe the dynamical nature of the molecular systems and thereby to provide a more realistic representation, which is often needed for understanding fundamental molecular properties. The method has proven to be remarkably successful in capturing molecular events and structural transitions highly relevant to the function and/or physicochemical properties of biomolecular systems. Herein, after a brief introduction to the method of MD, we use a number of membrane transport proteins studied in our laboratory as examples to showcase the scope and applicability of the method and its power in characterizing molecular motions of various magnitudes and time scales that are involved in the function of this important class of membrane proteins.

Living cells rely on the continuous exchange of diverse molecular species, e.g., nutrients, precursors, and reaction products, across the cellular membrane for their proper function.[1] Membrane transporters are proteins that serve as active gatekeepers closely regulating the traffic of these species across the membrane. Similar to channels, membrane transporters provide selective pathways for the permeation of their substrates across the membrane, but they are endowed with the additional ability to actively pump their substrates across the membrane, often against their electrochemical potential gradient. These specialized molecular devices provide the machinery to intimately couple active transport to various forms of cellular energy; light or the energy released by chemical reactions (mainly ATP hydrolysis) is utilized by primary active transporters, while the chemical gradients of other species are exploited by secondary active transporters.[1] The fundamental role of membrane transporters in diverse biological and physiological processes has rendered them as important drug targets, further stimulating widespread interest in mechanistic studies of these proteins at a molecular level.[2]

The widely accepted general mechanistic model for transporters, termed the alternating-access model,[3] proposes that the transporter protein switches substrate accessibility from one side of the membrane to the other, by undergoing structural transitions between outward-facing (OF) and inward-facing (IF) states, temporarily residing in several possible intermediate states. Crystal structures for several of these proposed states across various families of transporters have provided strong support for this mechanism.[4−7] However, the availability of structures for multiple functional and/or conformational states of the same protein, or even within the same family, is limited.[8−19] Even when multiple conformational states have been characterized for the same transporter protein, knowledge of the transitions between these states and the coupling of these transitions to the energy-providing mechanism is largely lacking for most membrane transporters. Hence, reconstructing the transport cycle and understanding the mechanism of energy coupling and substrate/ion transport remain challenging and rely on techniques that would yield a dynamical description of the process.

Computational techniques such as molecular modeling and simulation have been widely applied in the description of dynamics and mechanistic aspects of biomolecular function.[20] The overlap between computational and experimental observations has expanded considerably because of long time-scale simulations allowed by advanced, more efficient algorithms and the increased availability of faster and more powerful resources.[21,22] Computational methods are now routinely adopted for predictive or descriptive studies that parallel or complement experimental techniques and have allowed significant discoveries about how biomolecules work.

The application of computer simulation to studies of membrane transport proteins has a short history, because of the previously limited availability of structural information, as well as several technical challenges in computationally describing these complex systems.[20] These limitations have been addressed diligently over the past decade: advancements in crystallographic techniques have resulted in the determination of numerous crystal structures for several membrane transporters across diverse families,[4−7] and significantly longer simulations enjoying improved force fields, e.g., the recent revision of CHARMM providing more accurate representations of lipids and lipid–protein interactions,[23] have resulted in more realistic representations of physiologically relevant phenomena. These factors have together resulted in a quantum leap in the field of membrane protein simulations.[24−55]

Conformational changes are central to biomolecular function and particularly important to the mechanism of membrane transporters. The transport cycle of membrane transporters involves structural changes on a diverse scale, ranging from small (local) rearrangements such as salt bridge and hydrogen bond formation or breakage and localized gating motions to large (global) structural transitions such as helix or domain motions. The molecular events leading to these conformational changes and their effects on the overall mechanism are largely unknown and often difficult to characterize experimentally. Molecular dynamics (MD) simulation offers a method with sufficiently high temporal and spatial resolutions suitable for the characterization of conformational events ranging from local to global (Figure 1). Such structural transitions have been successfully captured in MD simulations of a large number of membrane transporters.[25−41,45−48,50,53,54,56−58]

Figure 1

Different scales at which structural and dynamical changes relevant to the function of membrane transporters occur. MD simulations can be used to study transporter dynamics at scales ranging from the global conformational and/or functional state of the protein (top) to changes in the secondary structure and topology (middle) to changes in local interactions between individual groups at the atomic level (bottom). The structure of LeuT is used to illustrate these levels. In the top panel, LeuT (white surface) is embedded in a lipid membrane (green), surrounded by water and ions. The crystal structure is in the OF state; the cavity through which the substrate binds is visible. The middle panel shows the LeuT secondary structure, where the 10 TM helices forming the unique LeuT fold are shown. In the bottom panel, specific salt bridge and hydrophobic interactions serve as extracellular gates in LeuT. Superimposed snapshots for the residues illustrate their dynamics observed during an MD simulation.

In this review, after providing a brief introduction to MD simulation, we describe the application of the method to membrane transporters, using examples from our own studies highlighting functionally relevant features captured by the simulations. We use five systems (Figure 2), viz., leucine transporter (LeuT), galactose transporter (vSGLT), glutamate transporter (GlT), glycerol 3-phosphate transporter (GlpT), and an ATP-binding cassette transporter (ABCT); this set samples mechanistically diverse classes of primary and secondary active transporters, the latter including symporters and an antiporter (Figure 2).

Figure 2

Five transporters discussed herein. The LeuT fold inverted repeat is highlighted in color in LeuT (green and blue) and vSGLT (green and pink). Bound substrate and ions are shown. The GlT trimer is shown, colored by the trimerization domain (orange) and transport domain (blue). The two helical bundles in GlpT are colored yellow and blue, respectively. For ABCT, the structure of the intact maltose transporter with different domains colored individually is shown. The approximate position of the embedding membrane is shown as two solid lines.

Molecular Dynamics Simulation

The molecular dynamics (MD) method and its applications to molecular systems are largely based on classical mechanics and statistical mechanical theories.[59] Intensive calculations are required in this method; hence, the efficiency of MD simulations is heavily dependent on algorithmic developments in mathematics and computer science as well as on advancements in computer hardware.[60]

In a typical atomistic MD simulation, interactions are calculated between atoms using a set of parameters that define a potential energy function, representing a “force field”. The derivation of force field parameters is typically a painstaking iterative process, in which initial parameters obtained from experimental data and quantum mechanical calculations[61] are optimized to reproduce structure and vibrational modes, as well as thermodynamic properties of the molecular systems of interest. Various force fields are available for biomolecular simulations, with minor differences in their potential energy functions, and corresponding differences in parameters.[62−65] A typical potential energy function for biomolecular simulations, U, includes terms that describe bonded (bonds, bond angles, dihedral angles, etc.) and nonbonded (van der Waals and electrostatic) interactions:Nonbonded interactions usually form the more expensive part of the calculations and often play a more important role than bonded interactions in describing interactions between different molecules, or even between different functional groups within the same molecule, e.g., side chain interaction within a protein.

During an MD simulation, at each time step, the total force F acting on each atom is calculated as the negative derivative of the potential energy U with respect to its coordinates r:Using these forces, the Newtonian equations of motion are then integrated[59,66] to determine changes in the position r and velocity v of individual atoms in time, t:where m is the mass and a is acceleration. The result of an MD simulation is, therefore, a trajectory file that includes an ensemble of configurations of the system over a period of time, typically ranging from tens to hundreds of nanoseconds for atomistic simulations of biomolecular systems. These trajectories contain information about motions and interactions that can occur in a system under given conditions and may also be used to calculate macroscopic properties for prediction of or comparison to experimental observables.

MD has the advantage of providing dynamical information at high spatial and temporal resolutions. It allows the researcher not only to monitor the natural motion of a molecular system in real time but also to test out various conditions, including those that are not accessible experimentally. While extremely powerful in generating and testing hypotheses, MD methods are severely limited in time scale because of the short time steps (on the order of femtoseconds for atomistic simulations) required for accurate integration of the trajectory and proper description of the molecular system. Many important biological processes occur on the microsecond, millisecond, or even longer time scales, while simulations for large biomolecular systems are currently limited to the nanosecond to microsecond time-scale range, hence restricting direct application of MD to the simulation of such processes completely. The time-scale limitation often translates into inadequate sampling of the configurational space and incomplete description of the entropy of the system. A second major shortcoming of classical MD simulations is the use of simplified potential energy functions to describe the molecular system that do not account well for electronic properties, particularly resulting in inadequate treatment of electronic polarization effects. To alleviate this problem, efforts aimed at developing “polarizable” force fields are under way.[67−70]

Despite these limitations, MD has been successfully employed in studying a wide range of biomolecular systems and phenomena, including membrane transport proteins and their mechanisms.[24,43,71] Furthermore, with continuous algorithmic improvement, the availability of faster hardware, and refinement of the available force fields addressing these limitations, the gap between simulations and experiments is rapidly closing, as evidenced by recent studies reporting simulations reaching time scales on the order of microseconds to submilliseconds.[72−75] We also note that often, long time-scale processes are long simply because they are composed of many “steps” or molecular events that each might occur in short time scales and therefore can be efficiently described using MD. Ideally, when structural descriptions of a sufficient number of intermediate states involved in a long process are available, MD can be used to reconstruct the dynamics of the entire process.

To study large transitions, as well as other high-barrier processes that are amenable to only prohibitively long equilibrium simulations, specialized simulation and sampling techniques that can be adopted to induce and study transitions within reachable simulation time scales have been developed. In techniques such as steered[76] and targeted MD[77] (SMD and TMD, respectively), external forces (or biasing potentials) are applied to one atom or a group of atoms, to induce specific molecular events and to drive the system from one state to another. Other enhanced sampling techniques that have been developed particularly with the purpose of calculating the free energy of processes such as substrate binding and release or conformational transitions include free energy perturbation (FEP),[78,79] thermodynamic integration (TI),[80,81] umbrella sampling (US),[82] and the more recently developed methods of adaptive biasing force (ABF)[83,84] and implicit ligand sampling (ILS).[85,86]

In a typical membrane protein simulation of the kind presented here, the researcher starts with an experimentally determined atomistic protein structure, to which water, lipids (membrane), and ions are added. Figure 1 shows a typical membrane protein simulation system, featuring LeuT embedded in a lipid bilayer.[41] NAMD[87] was adopted for the reported simulations. To mimic ambient experimental conditions, MD simulations are often conducted either under NPT [constant temperature (310 K) and constant pressure (1 atm), while keeping the total number of particles in the system constant] or under NVT conditions in which the volume of the system (instead of the pressure) is kept constant. Force field parameters are available for standard molecular systems, such as protein residues, nucleotides, and lipid molecules. Missing parameters and topology files for ligands are included either by adopting similar parameters from the available force field or by quantum mechanical calculations in a manner consistent with the employed force field. The simulations involve a brief period of initial membrane equilibration (typically 1–5 ns), wherein the lipid tails are allowed to equilibrate while the lipid headgroups and the protein are constrained to their initial positions. This is then followed by an unconstrained equilibration of the lipids and the protein for simulation times ranging from tens to hundreds of nanoseconds, dictated by the nature of the molecular problem and the available resources. As will be discussed, these time scales are appropriate for capturing various motions ranging from very fast ones, e.g., water and ion diffusion, to slower motions such as side chain reorientations and loop movements, and in some cases even larger conformational changes such as domain motions.

Applications of MD to Studying Membrane Transporters

Computational techniques often step in when a research problem becomes particularly refractory to experimental investigations. Prime examples of such research problems are the molecular details of the mechanisms involved in the function of membrane transporters. To understand the transport mechanism, a dynamical description of several states in the transport cycle, and the transition between them, is necessary. Despite recent progresses in the structural biology of membrane proteins, high-resolution structures of membrane transporters are still scarce, and multiple states are rarely captured for the same transporter. In some cases, the functional state represented by the determined structures is not straightforward to characterize. Furthermore, even in cases where multiple crystal structures exist for different states,[8−12,88] the dynamics of the transporter and molecular events involved in the mechanism, such as the sequence of binding, unbinding, and translocation events for the primary substrate and the cotransported ion(s), remain unresolved. Describing such mechanistic aspects requires methodologies offering a dynamical treatment of the protein of interest. Recent applications of experimental methods such as electron paramagnetic resonance (EPR) and single-molecule Förster resonance energy transfer (FRET) have shown promise in characterizing the conformational dynamics of membrane transporters.[89−92] Further application of these techniques to multiple transporter families is expected to greatly enhance the current understanding of the mechanism of transport. Meanwhile, MD simulations have proven to be a powerful tool, complementary to experimental approaches, for characterizing transporter dynamics.

In the context of membrane transporters, some of the problems that can be investigated by MD simulations are the study of atomic-scale details of interactions of the protein, substrate, ions, and water; substrate–ion coupling and cotransport; the nature of permeation and binding–unbinding pathways; the effect of titration states of protein residues or the substrate on the transport mechanism; and substrate-induced protein conformational changes. Moreover, in conjunction with other modeling techniques, MD can be employed to provide structural models for missing states in the transport cycle and to provide a possible mechanism for structural transitions between major states.[25,41,93,94] Finally, MD simulations provide an extremely powerful method for characterizing the energetics associated with various processes and transitions in biomolecular systems. Such calculations are among the most expensive computations but have been successfully employed, e.g., to compare the affinity of binding of ions to various binding sites in transporters[26,95,96] and to calculate the free energy profile of substrate translocation along its binding pathway.[97,98]

MD has been most effective in capturing small-scale structural changes in membrane transporters, such as gating motions caused by side chain flipping, loop motions, and the movement of small molecular species, e.g., water, ions, and, in some cases, the primary substrate itself. The study of larger-scale changes in transporters using MD is restricted by the time scales that can be reached using conventional MD methods. This restriction, coupled with the stochastic nature of conformational changes, considerably lowers the probability of capturing global structural changes in MD simulations. While in rare cases large-scale changes have been observed within a few hundreds of nanoseconds in equilibrium MD simulations of membrane transporters,[40,93] such transitions may be studied by employing biased simulations to accelerate the process, e.g., by performing SMD or TMD.[25,31,43,97,99]

The following sections will describe the application of MD and other modeling techniques to the mechanistic investigation of these five transporter systems, viz., LeuT, vSGLT, GlT, GlpT (all secondary active transporters), and an ABCT (a primary active transporter). LeuT, vSGLT, and GlT perform Na+-coupled symport, while GlpT utilizes a substrate gradient to perform antiport. ABCTs are primary active transporters that utilize ATP binding and hydrolysis as a source of energy for substrate transport. The major focus of each section will be the design of the computational approach and its application, as well as the nature of the lessons learned from the application of the methods. We have therefore refrained from providing a detailed description of specifics of each molecular system and the results. The reader is encouraged to refer to the original articles for more details about each system. The spectrum of the conformational changes characterized in these studies is presented as two classes: smaller-scale events that are often efficiently sampled by MD simulations and large-scale motions that rely on extended equilibrium simulations and/or biased simulations. In the accompanying figures, these two types of molecular motions have been distinguished by using green and blue frames, respectively.

Modeling an Unknown State for Neurotransmitter Transporters

The leucine transporter (LeuT) is an amino acid transporter, and a bacterial homologue of the neurotransmitter sodium symporter (NSS) family that includes transporters responsible for clearing the synaptic cleft of neurotransmitters, e.g., serotonin, dopamine, norepinephrine, and GABA. The NSS family thus includes important drug targets for the treatment of neurological conditions such as depression, anxiety, and drug addiction.[100−102] LeuT provides the only crystal structure known for this family and also represents the first example of a unique fold, termed the LeuT fold,[18] characterized by a 5+5 inverted-repeat arrangement of TM helices and the presence of two broken helices in the center. The LeuT fold was subsequently observed in several other families of transporters, suggesting a possibly common mechanism of transport despite functional diversity.[6,7,103] LeuT thus serves as an important model for the study of the dynamics and mechanism not only for NSS members but also for the whole family of LeuT fold transporters.

Despite having access to structural information for only the OF state[18,104−110] until recently,[9] several studies focused on deciphering the transport mechanism of LeuT have been reported over the years. The field has particularly benefited from several collaborative studies between experimental and computational researchers jointly reporting models of the IF state[25,111] and the dynamics of LeuT relevant to its transport function.[89,90,112,113]

LeuT function is known to be Na+-dependent, but the structural and dynamical roles of Na+ ions in the transport function have emerged only recently. Several simulation studies of LeuT have addressed the dynamics of the two Na+ ions bound to LeuT in most of its reported structures.[114] Earlier MD studies coupled with free energy calculations have revealed that the bound Na+ ions stabilize the binding pocket and that both ions were required for effective substrate coupling.[115] Structural and energetic factors determining ion selectivity in LeuT have also been examined.[26,116] An interesting observation noted in MD studies of LeuT, its homologue SERT, and other LeuT fold transporters, vSGLT and Mhp1, has been the spontaneous release of a bound Na+ from the IF state.[8,25,35,41,93,117] This event occurs while the substrate is still bound, suggesting the possible order of release of ions and substrate to the intracellular side. Furthermore, in simulations of substrate-free LeuT with bound Na+, and in corresponding EPR studies, the presence of Na+ was observed to increase the accessibility of the binding site that could aid substrate binding, thus suggesting the possible order of binding of ions and substrate.[90]

The earliest MD studies of LeuT used the OF crystal structure to describe the putative substrate binding–unbinding process and the transport mechanism.[25,31] SMD simulations in which forces were applied on the substrate bound to LeuT were performed to induce its release to the extracellular solution.[25,31] During the induced translocation of the substrate through the extracellular vestibule, major barriers to unbinding were offered by interaction with a bound Na+ ion, aromatic residues (Y108 and F253) forming an occlusion over the substrate, and a salt bridge (R30 and D404). Rotation of the side chains of Y108 and F253 was found to allow opening of the binding site when the substrate exited. The substrate showed a common unbinding pathway in a majority of the simulations, suggesting the functional relevance of this pathway.[31] SMD studies of LeuT also revealed for the first time the location of a putative secondary substrate binding site, which was later verified experimentally and became a topic for much debate over its significance to the transport mechanism.[18,25,105−107,112,118,119]

Over the past few years, several crystal structures of LeuT have been reported, all representing only the OF state,[18,104−110] until a recent report in which the IF state structure was captured.[9] While the lack of structures of IF and intermediate states continued for several years, multiple attempts were made to model these alternate states of LeuT.[25,41,111] The absence of crystal structures for homologous proteins presented a challenge as standard modeling techniques could not be adopted to generate model structures of the unknown states. In one of the first computational studies aimed at modeling the unknown LeuT IF state, a novel approach in which the structure of one-half of the inverted repeat was switched with the other half, resulting in a structure that was open to the intracellular side, was adopted.[111] Because the two inverted repeats lacked significant homology, accurate positioning of residues could not be expected; however, the TM helix positions provided mechanistic insight. Most significantly, the interconversion of the OF and IF states was attributed mostly to the rocking of the TM1–TM2/TM6–TM7 bundle as a rigid body with respect to the rest of the protein. In another approach, the substrate was pulled through LeuT toward the intracellular side using SMD simulations, to induce inward opening.[25] Because this study used the LeuT crystal structure as the starting conformation, it could record residue–residue interaction changes and provided information about local interaction network rearrangements that are possibly involved in the formation of the IF state. In a later study using extended MD simulations, a comparison of the LeuT interaction network with that in the IF state of ApcT, a related LeuT fold transporter, provided details about the common mechanistic aspects of IF state formation in the LeuT fold.[25,117] However, the SMD-based modeling method was limited in capturing larger global conformational changes that might arise during the transport process.

In a third approach, a modeling methodology in which the aim was to generate the IF state with much of the structural information of the OF state intact was devised, while also attempting to capture local and global structural changes occurring during the formation of the alternate state. The methodology combined computational techniques of structure-based alignment, comparative modeling, targeted MD, and equilibrium MD.[41] The structure of vSGLT,[120] a LeuT-fold transporter in the IF state, was used as the template for building the IF state model. Structure-based alignment was performed to align the structurally similar portions (5+5 TM repeats) between LeuT and vSGLT, and the resulting alignment was used to generate a preliminary model of IF state LeuT. Rather than this homology model being used directly as an IF structure for LeuT, TMD was used to gradually deform key helices in the OF crystal structure toward the conformation in the IF homology model. The rest of the protein and side chains were allowed to undergo free dynamics in their natural environment of water, ions, and membrane, hence incorporating the effect of the environment on the protein as the transition was induced[41] and ensuring a better representation of the complete protein structure in the resulting model. Equilibrium MD was then performed to obtain a relaxed model of the IF state structure.

During the TMD simulations in which the OF to IF transition was induced, several local interactions were broken in the intracellular half, while new ones formed in the extracellular half. In particular, the rearrangements of hydrogen bonds, salt bridges, and hydrophobic interactions were monitored through the simulations (Figure 3A). Conformational hot spots were thus pinpointed, revealing the structural elements of the transporter that contributed to the transition (Figure 3A).[41]

Figure 3

Structural changes of LeuT during the OF to IF transition. (A) Salt bridge rearrangements captured in the extracellular (EC) and intracellular (IC) halves of LeuT as the structure transits from the OF to IF state. Side chain positions before (light) and after (dark) the transition are shown. Distances between salt bridging residues show salt bridge formation as closure occurs in the EC half and breakage as opening occurs in the IC half. (B) OF to IF state transition in LeuT. Superimposed structures before (transparent) and after (colored) the transition viewed along the plane of the membrane (left) and perpendicular to the membrane: the extracellular half (right, top) and intracellular half (right, bottom). A contact map of interactions broken (black) and newly formed (red) during the transition, pinpointing conformational hot spots (green ovals).

The TMD simulations designed for IF model generation also allowed the monitoring of global conformational changes that transform the structure from the OF state to the IF state (Figure 3B).[41] It should be noted, though, that the transition pathway observed in such biased simulations does not necessarily represent the biologically relevant transition pathway, and further examination of the pathway using parallel computational or experimental techniques is recommended. Experimental mutagenesis and accessibility studies using MTS reagents have revealed that the cytoplasmic lumen in SERT, a mammalian homologue of LeuT, is lined by residues from TM1, TM5, TM6, and TM8.[111,121] The TMD-based LeuT IF model successfully captures these observations. Recent reports of LeuT dynamics using FRET and EPR have identified important dynamic components, including the role of the outward swinging of TM1 in the formation of the IF state.[89,90] The dynamics of LeuT monitored during simulations agreed with these observations. These simulations also yielded mechanistic details that have not been characterized before: in the IF state, Na+ unbinding was observed to occur before the substrate, suggesting the sequence of events in the transport cycle; the TM1–TM2/TM6–TM7 bundle was found to be internally flexible, i.e., did not move as a rigid body as suggested in an earlier model;[111] and a possible role for TM2 and TM7 as facilitators was also predicted. Detailed contact maps (Figure 3B), indicating the changes in interactions from the OF to IF state, were generated using the simulation results, providing an informed starting point for further mutagenesis studies designed to study the LeuT transport cycle.

The much-awaited IF state crystal structure of LeuT was reported recently at 3.2 Å resolution.[9] The structure broadly confirms some predictions about the overall structure of the IF state from the three modeling attempts described above,[25,43,111] such as the helices lining the intracellular opening and the swinging out of TM1 as an important feature of inward opening. However, the extent of TM1 swinging and the limited role of other TM helices in comparison to that of TM1 evidenced by this structure had not been addressed by these studies. The first methodology in which the two halves of the inverted repeat were switched[111] is computationally inexpensive and yields a coarse model that can be used to study the secondary and tertiary structural arrangement. The other two methods, based on SMD[25] and TMD,[43] respectively, are time-consuming but yield a much higher level of information in terms of possible residue–residue interaction changes, dynamics of the ions and substrate, changes in solvent accessibility, etc. Structural superimposition of the freely available model generated by the TMD-based method[43] on the crystal structure resulted in a backbone root-mean-square deviation (rmsd) of 3.4 Å for all residues, which is reduced to 2.8 Å upon exclusion of TM1. Considering that in nanosecond-scale simulations of membrane proteins with low-resolution crystal structures, the rmsd often reaches 3.0 Å, the model performs fairly well in this comparison.

The LeuT IF state modeling efforts thus provide an example in which modeling and simulation techniques were used effectively to study local structural rearrangements (Figure 3A), as well as global conformational transitions (Figure 3B), in membrane transporters. Also, the power of a combination of approaches such as homology modeling and biased MD with conventional MD, when any single computational technique is insufficient, is demonstrated.

Sequence of Unbinding Events in Solute Sodium Symporters

The solute sodium symporter (SSS) family consists of secondary active transporters that couple Na+ symport to the transport of a wide range of solutes, including sugars, amino acids, vitamins, and inorganic ions.[122] Members of the SSS family play crucial roles in human health, and their malfunction is associated with various metabolic disorders.[123] vSGLT is the first SSS transporter for which high-resolution structures have been determined.[37,120] Although it adopts the LeuT fold, it is functionally divergent from LeuT.

MD studies of vSGLT have provided much insight into key aspects of the transport mechanism, i.e., its dynamics and transport mechanism,[35−37,50] most prominently characterization of the functional state of the crystallographically captured structure, the functional importance of the structural element missing in the crystal structure,[124] and permeation of water through vSGLT.[125,126]

In the first crystal structure of vSGLT, which was reported to be a substrate-bound IF state,[120] the position of the Na+ ion could not be directly resolved, but on the basis of the structural homology to LeuT for which the position of the ions had been unequivocally determined,[18] a Na+ ion was also modeled in the reported Protein Data Bank file for vSGLT.[120] In the crystal structure, the substrate, galactose, is bound approximately halfway across the membrane and proposed to be occluded from the cytoplasmic side by hydrophobic residues, especially Y263 stacking with the pyronose ring of the substrate (Figure 4A).[120] Simulation studies, however, strongly suggested that the state captured in the crystal structure was indeed representing a Na+-free (Na+-releasing) structure,[35,36] likely representing an IF-open state rather than an IF-occluded state.[50]

Figure 4

Cytoplasmic release of the substrate and the cotransported ion in a solute sodium symporter. (A) Overview of the substrate release trajectory shown during equilibrium MD (pink) and SMD (red). Residues lining the substrate pathway are shown as sticks. The substrate (galactose) and the suggested gate residue, Y263, are shown in van der Waals form. (B) Time evolution of the z coordinate of the Na+ ion during six independent simulations highlighting its spontaneous unbinding from the putative site in the crystal structure. The blue bar highlights the region in the vicinity of D189. (C) Position of the substrate during a 200 ns equilibrium simulation shown using the maximal and minimal z coordinates (blue solid lines) and the geometrical center (black solid line) of the substrate. For comparison, the z position of the geometrical center of the ring of Y263 (red solid line) is also shown. Two arrows highlight the two observed substrate unbinding events from the binding site.

In a series of equilibrium simulations (10 ns each) of vSGLT, starting from the crystal structure described above,[120] surprisingly, rapid unbinding of Na+ from its reported binding site was consistently observed in multiple independent runs (Figure 4B).[35] This functionally relevant event was also observed subsequently in simulation studies reported by other laboratories.[36,37] On the basis of these simulations, a low-affinity Na+ binding site was proposed for the IF state of vSGLT, an attribute that is strongly corroborated by structural comparison of the site with that in the OF state of LeuT.[35] Given that in these simulations the protein did not undergo a major structural deviation from the reported crystal structure (rmsd of ≤2.5 Å), the observed spontaneous Na+ release suggests that the crystal structure of vSGLT is in a Na+-releasing state.[35]

During its release, the ion hovers around an aspartic acid, D189, ∼4.5 Å from the proposed Na+ binding site in the crystal structure. The interaction between the Na+ ion and D189 was consistently observed in the MD simulations.[35−37] The functional importance of D189 is supported by experimental studies of homologous proteins, hSGLT1 and Na+/proline transporter PutP, where mutation of the corresponding residue results in complex effects ranging from a significant decrease in the level of Na+-coupled transport[127] to altered Na+ versus H+ selectivity,[128] supporting the notion that this highly conserved residue might play a role in Na+ unbinding in SSS transporters.

Subsequent extended equilibrium simulation (200 ns) revealed also spontaneous unbinding of the substrate from the conformational state captured in the crystal structure (Figure 4C).[50] In contrast to the previously proposed gating role of a tyrosine (Y263) side chain,[36,37] the unbinding mechanism captured in our equilibrium simulation does not rely on the displacement and/or rotation of this side chain. Rather, during its unbinding, the substrate adopts a curved pathway going around the side chain of the proposed gate, thereby avoiding the need for its rotation. The revealed pathway, which can be easily obfuscated from plain visual examination of the static crystal structure, could be revealed only through the extended simulations in which natural thermal fluctuation of the substrate within its binding pocket, assisted by stepwise hydration of the binding site, resulted in the spontaneous movement of the substrate along this hidden pathway.[50]

During the spontaneous exit of the substrate from the binding site, H-bonds with T431 and N142 appear to facilitate substrate unbinding.[50] As support for their involvement in the exit pathway, it has been shown that the mutation of the residue corresponding to N142 in the homologous human protein SGLT1 impairs transport.[129] The residue corresponding to T431 in vSGLT (T460 in human SGLT1) was even suggested to be a substrate-binding residue in SGLT1 because its mutation to cysteine altered sugar selectivity and decreased the affinity for glucose.[130] Interestingly, neither N142 nor T431 appears to contribute directly to the substrate binding site in the crystal structure,[120] and only the unbinding mechanism and pathway characterized by the simulations could provide a molecular explanation for their importance for proper transport function.[50]

Although the substrate exhibits full unbinding from its crystallographically determined binding pocket on several occasions during the equilibrium simulation, it does not completely leave the protein’s lumen during the simulated time scale, and thus, the unbinding event described above was only partial. To probe the remainder of the exit pathway from the protein lumen into the cytoplasmic milieu and to examine the existence of other barriers along the pathway, SMD simulations in which the substrate was slowly pulled toward the cytoplasmic solution were employed (Figure 4A).[50] These biased simulations were repeated each time starting from a different snapshot taken from the equilibrium trajectory representing various configurations of the substrate sampled in the absence of any forces. The results confirmed that once the substrate avoids the major obstacle of Y263, an event observed during the equilibrium simulation, it can smoothly diffuse along the rest of the exit pathway toward the cytoplasm without facing any major barriers and, importantly, without requiring any conformational change in the transporter protein. The observed molecular events clearly indicate that no gating motion is required for the release of the substrate from the crystallographically captured structure of the IF state of vSGLT.[120] In other words, this conformation likely represents an open, rather than occluded, state of the transporter.[50] Several residues, e.g., E68, N142, T431, and N267, facilitate the substrate’s initial displacement from the binding site, while along the remainder of the exit pathway, the translocation of the substrate is facilitated by the sequential formation and breakage of H-bonds between the substrate and a series of conserved, polar residues lining the pathway.[50]

Computational studies of vSGLT[35−37,50,124−126] provide a typical example of using MD simulations to access the dynamics of small molecules such as ions, water, and the substrate and thereby to characterize the functional state of crystallographically determined structures of membrane transporters, a task that is often not straightforward.

Gating and Substrate–Ion Coupling in the Glutamate Transporter

GlTs [also termed excitatory amino acid transporters (EAATs)] belong to solute carrier family 1 (SLC1) of neurotransmitter transporters, which catalyze translocation of neurotransmitters from the synapse to the cell.[131] The malfunction of GlT in humans is associated with numerous diseases and pathological states, including neurodegenerative disorders, epilepsy, schizophrenia, and traumatic brain injury.[132−134] Driven by pre-established electrochemical ionic gradients, GlT catalyzes the “uphill” translocation of its substrate. During each transport cycle, one substrate is translocated along with three Na+ ions. The crystal structures of a bacterial homologue of GlT (Gltph) in both the IF and OF states[11−14] have provided a structural framework for understanding the transport mechanism of all GlTs.

These structures have revealed a trimeric architecture for GlT along with the positions of the substrate and two of the Na+ binding sites (termed Na1 and Na2 sites) in each monomer. Although the structures do not provide any information about the third Na+ binding site (Na3), mutagenesis studies[135] of a mammalian GlT have provided hints about the putative Na3 site. It has been shown, for instance, that an acidic residue corresponding to Asp312 in Gltph, which is not coordinated to any ion in the crystal structure, is involved in binding to one of the Na+ ions during the transport cycle, suggesting that this residue is likely involved in the Na3 site. Various experimental studies[135−137] have indicated that the binding of the substrate or Na+ ions induces conformational changes in GlT that facilitate the movement of these species inside the transporter. These studies clearly suggest a coupling between the binding and translocation of these species during the transport cycle. However, their limited resolution did not allow one to draw specific conclusions about the nature of such conformational changes and substrate–Na+ coupling.

To investigate the effect of binding of the ions and the substrate on protein dynamics, one would ideally design and perform equilibrium MD simulations of the protein in the presence of these species initially placed in the bulk and monitor their movement over time as they diffuse, find, and bind to their respective binding sites. While such calculations can, in principle, be done (e.g., see refs (40), (55), and (138−140)), and one might even increase the likelihood of capturing successful events by including a large number of ligands in the simulations (flooding),[29,86,141] the time required for free diffusion of the substrate (or ions) from the bulk into the lumen of the protein and then into their respective binding sites is usually prohibitively long, especially when the binding site of interest is buried inside the protein. Because it was not expected that the complete binding of the substrate and Na+ ions could be fully captured by equilibrium MD simulations, a different strategy was taken for the case of Gltph. Starting from the OF crystal structure,[14] we performed multiple simulations, in which the substrate and/or different Na+ ions were either kept or removed at the beginning of each simulation, thereby generating a set of simulations that included all possible combinations of the protein, the substrate, and the ions.[32,39] The goal was to study the dynamics of the protein in these different “bound” states and to use the resulting dynamical information to draw conclusions about the nature of the coupling between the various elements and, we hope, about the sequence of binding events. These simulations were performed on membrane-embedded trimeric models, which allowed us to simulate three different bound states within each simulation systems, as the individual monomers are believed to function independently (no cooperativity).

Monitoring structural fluctuations of Gltph in these simulations revealed the gating dynamics of this transporter.[14,27,32] Comparison of the dynamics of the substrate-bound and substrate-free simulations showed a clear contrast in the dynamics of helical hairpin HP2 in these states (Figure 5). In all the simulations performed in the presence of the substrate, HP2 invariably stays closed, while removing the substrate results in a large opening motion of HP2 and complete exposure of the substrate to the extracellular solution (Figure 5A,B).[32] These results suggest that HP2 plays the role of the extracellular gate and, more importantly, that its opening and closure are controlled by substrate binding.[32] A gatinglike motion of HP2, which has also been observed in an MD study by another laboratory,[27] is supported by the crystal structure of Gltph in the presence of an inhibitor.[14]

Figure 5

Dynamics of the extracellular gate and substrate–ion coupling in the glutamate transporter. (A) HP2 loop motions responsible for extracellular gating. Superimposed snapshots of HP2 (gray) show opening motions after the removal of the substrate. Yellow and pink show the closed conformation of HP2 and its open conformation in one of the last simulation frames, respectively. (B) Time evolution of the rmsds of helical hairpins HP1 and HP2 in the presence (top) and absence (bottom) of the substrate. (C) Substrate-induced formation of the Na2 site (marked with a red oval). The dipole moments (arrows) of two short helices, from TM7 and HP2, which are misaligned in the substrate-free state, converge at a point upon substrate binding, creating a highly negative electrostatic potential for binding of the Na2 ion. (D) Putative site for Na3. The residues involved in direct coordination of the ion in this site are labeled.

The simulations also revealed that substrate binding results in the formation of one of the Na+ binding sites (Na2 site) at a position between two short helices.[32] In the apo state, the dipole moments of these helices are not aligned. When the substrate binds, these two opposing helices align such that their dipole moments converge on a common spot, resulting in the formation of the Na2 site (Figure 5C). These results suggest that the formation of the Na2 site is likely an event induced by the binding of the substrate and that substrate binding precedes that of the Na2 ion. These aspects are both in agreement with experimental current measurements, indicating that substrate binding allows the binding of one of the Na+ ions.[142,143] In addition, the functional role of Na1 in Gltph has been described in another simulation study[144] employing metadynamics, a nonequilibrium simulation technique used to reach larger-scale conformational changes.[145−149] Metadynamics was used to describe the transition of the OF state between the open and occluded state.[144] Although using a single Gltph monomer immersed in water[144] represents a drawback in these simulations, the results indicated that binding of Na1 stabilizes the open conformation, which is consistent with its facilitatory role in substrate binding.[32,135]

To probe potential binding sites for the third Na+ ion (Na3), for which no binding site has been identified in the crystal structures, we took advantage of the fast dynamics of water molecules because they can quickly visit sufficiently polar regions within the protein matrix. We focused our analysis on water molecules that hydrate the region around Asp312 during the simulations, as this residue has been implicated in ion binding based on mutagenesis studies[135] showing that neutralization of this residue inhibits binding of Na+ to the substrate-free form and thereby prevents the cycling of the glutamate transporter. A series of MD simulations were performed, in which individual water molecules visiting this region were selected one by one and “mutated” to a Na+ ion. The dynamics of the inserted Na+ ion and its interaction with the surrounding site were monitored,[39] to identify a consensus site to which Na+ ions frequently bind. A putative Na3 site was thus identified, in which the Na+ ion is coordinated by Asp312, Asn310, Thr92, and one water molecule (Figure 5D).[39] The putative Na3 site suggested by the simulations is strongly supported by recent mutagenesis experiments with EACC1[150] showing that the T101A mutation (corresponding to Thr92 in Gltph) dramatically decreases the apparent affinity of Na+ for empty EAAC1, indicating the direct involvement of this conserved residue in Na+ binding in GlTs. In a different computational study,[96] a different binding site for Na3 was proposed on the basis of grand canonical Monte Carlo and MD simulations, which is ∼8 Å from our proposed Na3 binding site and is coordinated by Thr314, Ala353, and Asn401 as well as the bound substrate. This binding site was also supported by mutagenesis studies showing that mutations of Thr314A and Asn401A abolish Na+-driven glutamate uptake.[96] Although the two binding sites are different, their relevance to the transport process, as indicated by experiments, might indicate that they may become important at different stages of the transport cycle.

The described simulations for a bacterial GlT captured some of the functionally relevant features for this transporter, including the nature and motion of the extracellular gate at an atomic level, how its closing and opening are controlled through substrate binding, and the binding site for Na3. The dynamics of the IF state of Gltph has also been studied with MD simulations,[46] suggesting the involvement of both hairpins HP1 and HP2 in intracellular gating. Moreover, the global conformational changes during the transition between the OF and IF states have been explored using the anisotropic network model,[151] revealing two major types of large-scale motions: an asymmetric stretching and contraction and a symmetric opening and closing of the three subunits. These theoretical studies, together with the experimentally characterized mechanistic aspects,[14,135−137] provide a more complete picture of the transport mechanism in GlTs.

Substrate-Induced Rocker Motion in an Antiporter

GlpT belongs to the major facilitator superfamily (MFS),[152,153] the largest superfamily of secondary active transporters. Found in all kingdoms of life, MFS transporters include several medically and pharmacologically important proteins, e.g., efflux pumps conferring resistance to antibiotics in bacteria and to chemotherapeutics in cancer cells.[152] GlpT, which couples the import of glycerol 3-phosphate (G3P) across the cytoplasmic membrane to the export of inorganic phosphate (Pi), was among the first structurally characterized members of the MFS.[153−158] Thus, the crystal structure of GlpT has been also used as a structural and functional model for other MFS proteins.[152,153]

MFS transporters typically consist of a single polypeptide chain folded into a pseudosymmetric pair of helical bundles each composed of six transmembrane α-helices.[153,159,160] The crystal structures available for MFS transporters[153−158] support a “rocker-switch”-like motion of the two bundles during the transport cycle.[153,161−163] Following the proposal[164] that MFS transporters share the inverted repeat topology observed in other transporter families,[165,166] a structural model of lactose permease (LacY) from MFS in the OF conformation was obtained by swapping the conformations of inverted-topology repeats found in the protein, which conforms to the alternating access model both theoretically and experimentally.[164]

GlpT has been crystallized without a bound substrate, i.e., in its apo state,[153] and only an approximate putative binding site for the substrate could be inferred from the structure. As substrate binding is believed to be coupled to protein conformational changes relevant to the transport cycle, describing the binding mode of the substrate in a detailed manner is one of the first steps in a mechanistic characterization of GlpT.

Docking is one of the most popular computational approaches used for identifying potential binding sites and favorable ligand–substrate poses. In docking, the favorable locations, orientations, and conformations of the ligand within the protein of interest are probed through efficient search algorithms and ranked on the basis of scoring functions.[167,168] Despite computational efficiency, most docking procedures at best only partially take into account the intrinsic flexibility of the protein, an important component that we have only recently been able to start to take into consideration.[169,170] Furthermore, docking cannot provide information about the nature of the binding pathway and mechanism. MD simulations, on the other hand, allow for direct observation of the binding process and a less biased description of binding sites and modes. They are, however, computationally much more expensive, with a nontrivial chance of not yielding any answer to the problem, as diffusion and binding of the substrate to its respective binding site might take a long time. Given the highly charged nature of the substrate(s) in GlpT, and the expectation that strong electrostatic forces might be involved, we decided to test whether equilibrium MD simulations might be able to capture the binding mode and site relying only on unbiased diffusion of the substrates from their initial position at the mouth of the lumen to the binding site.[33,40]

Another important aspect of these simulations, which is often important for other biomolecular problems, is the choice of the protonation states of titratable groups. Conventional MD simulations cannot accommodate titration state fluctuations in response to changes in the local environment. The majority of the reported MD studies, including our own, therefore, have considered only a “fixed” set of protonation states for the titratable groups of the molecular systems of interest. To set these protonation states during the initial setting of the system, a variety of algorithms (available in different programs and through various web interfaces) have been developed that allow the user to quickly estimate pKa values of residues in the presence of the rest of the molecular system.[171−179] Once set, the protonation states of the residues will not change over the course of the MD simulations, even though significant structural changes that might alter the local environment of these residues might indeed occur, in particular in extended simulations. To alleviate this problem to some degree, one might reassess the protonation states at certain intervals and if needed readjust them before continuing the simulations. A more systematic treatment is provided by the so-called “constant-pH” simulation protocols that allow a more dynamical, on the fly assignment and adjustment of the titration states of chemical groups.[180−184] These simulations are, however, much more expensive and do not yet fully cover the wide range of titratable groups available in proteins, lipids, etc.

For GlpT, one of the main groups that was expected to be heavily affected by the problem of fixed protonation described above is the substrate itself. Having in mind the high efficiency of running multiple parallel simulations when one uses a large number of processors, we opted to treat each physiologically relevant titration state of the substrates of interest in separate simulations.[33,40] The second pKa value for the substrates of GlpT (Pi and G3P) is very close to neutral pH; that is, monovalent and divalent forms of both substrates have approximately equal concentrations under physiological conditions. To account for this, multiple simulations each with a different titration state of each substrate (Pi–, Pi2–, G3P–, and G3P2–) were performed.[33,40] These simulations consistently revealed both local and global conformational changes that accompany substrate binding, thus characterizing substrate binding and protein dynamics simultaneously (Figure 6).

Figure 6

Spontaneous substrate binding and substrate-induced conformational changes in GlpT. (A) Trajectory of Pi binding (left) starting from its initial position (van der Waals representation, red) to its final binding site (blue). K80 acts as a “fishing hook”, catching the substrate at the mouth of the lumen and escorting it to the binding site at the apex of the lumen where it mainly interacts with R45 and H165. Distances between Pi and key residues in the binding site (right). No direct contact between Pi and R269 is observed during the simulations. (B) Substrate-induced helical motion in GlpT. H5 and H11 become straight upon substrate binding on the cytoplasmic side, resulting in partial occlusion of the transporter in this region. Three superimposed conformations of H5 and H11 (left): initial (t = 0; black) and final (t = 50 ns; red) along with an intermediate snapshot (t = 25 ns; blue). Distances between the H5 and H11 helices (right) measured by Cα atoms of representative residues located on the same x–y plane. All substrate binding simulations result in closure of the cytoplasmic side, while apo system simulations do not. Different substrate binding simulations are colored differently.

Experimental mutagenesis of highly conserved residues in the putative binding site (R45, R269, H165, and K80) results in the loss of binding and transport.[56] Moreover, unlike all other arginine residues, which could be mutated to lysines without impairing function, residues corresponding to R45 and R269 in a homologous protein (UhpT) were shown to be vital for the function.[153] The spontaneous substrate binding observed consistently in all the simulations not only supports the binding site proposed on the basis of the crystal structure but also demonstrates the dynamics and distinct role of each residue.[33,40] The simulations revealed that K80 functions as a “hook” in recruiting and escorting substrates from the mouth to the apex of the lumen, aspects consistently observed in all the simulations (Figure 6A). The phosphate moiety rapidly associates with R45 in the apex whose guanidinium group acts like a “fork” stabilizing the substrate in its binding site. Moreover, the hydroxyl groups of three conserved tyrosine residues (Y38, Y42, and Y76) residing below the guanidinium group of R45 form a “cage”-like binding pocket providing hydrogen bonding. The role of these conserved tyrosines was previously suspected to be only maintaining the basicity of the binding site.[56,152] The simulations, however, characterized this “tyrosine cage” as a part of the binding site. R269 was also suspected to contribute to substrate binding in a manner similar to that of R45, because of their symmetric arrangement in the crystal structure. However, R269 did not directly bind the substrate in any of the simulations.[33,40] Given the experimentally established importance of R269 for function,[56] the absence of its direct interaction with the substrate immediately after substrate binding as observed in the simulations suggests that it might play a role in later stages of the transport cycle.[40] Furthermore, the substrate failed to diffuse into the binding site upon in silico neutralization of R269, suggesting that the positive charge of this residue makes an important contribution to the attractive positive electrostatic potential of the lumen.[40]

Several experimentally testable hypotheses about the substrate binding site emerged from these simulations. Two distinct components of substrate–protein interactions could be identified: phosphate–GlpT interaction identified in both Pi2– and G3P2– simulations and glycerol–GlpT interactions identified only in the G3P2– simulations.[33] Binding affinities measured experimentally for various mutants clearly verified the predictions made on the basis of the simulations; while mutation of the residues proposed to bind the phosphate moiety disrupts the binding of both Pi and G3P, mutation of glycerol-interacting residues disrupts only G3P binding,[33] with no effect on Pi binding. Therefore, the first component of the binding site accounts for general recognition of phosphate-containing species, while the second component confers higher binding affinity to G3P.[33]

In addition to the local conformational changes during substrate recruitment, several salt bridge rearrangements also occur as a result of substrate binding. Salt bridges play important roles in transporter conformational changes likely by stabilizing different functional states. The closed periplasmic side of GlpT is kept together by a salt bridge network formed by K46 on the N-terminal half and D274 and E299 on the C-terminal half. These residues are necessary for transport but do not contribute directly to substrate binding.[56] The simulations show that substrate binding disrupts this salt bridge network, primarily by manipulating K46. Perturbation of the salt bridges upon substrate binding has been also reported for the mitochondrial ADP/ATP carrier (AAC),[138,139] suggesting a general mechanistic element for how substrate binding, particularly for charged substrates, might facilitate the progression of larger-scale conformational changes.

Simulations of GlpT provide one of the few examples in which functionally relevant large-scale motions were also captured upon substrate binding.[40] Two initially bent helices (H5 and H11) become straight upon substrate binding on their cytoplasmic sides, partially closing the mouth of the lumen (Figure 6A).[40] Interestingly, these helices appear to be inherently more flexible than the other helices in GlpT when simulated individually.[28] These helical conformational changes were observed only in the presence of the substrate, i.e., not in simulations, indicating their substrate-dependent nature.[40] The resulting structure is partially closed to both sides of the membrane[40,153] and, therefore, reminiscent of a “substrate-occluded” state inferred from the “alternating-access mechanism”.[3]

Simulations described in this section elucidated several functionally relevant molecular events and provided a dynamical view for an MFS transporter. These events include the coupling between substrate binding and the onset of protein conformational changes along the required alternating-access mechanism. These observations detailing the molecular interactions have been strongly supported by the experimental measurements designed to test and verify them.

Conformational Response of ABC Transporters to ATP Hydrolysis

ABC transporters make up a superfamily of primary active transporters found in almost all organisms and responsible for essential processes such as acquisition of nutrients, removal of hazardous molecules, and homeostasis of a wide variety of compounds.[185] Both normal activity and malfunction of different ABC transporters have been directly implicated in many diseases, most prominently in cystic fibrosis,[186] and in the development of multidrug resistance in cancer cells and pathogenic organisms.[187]

Most functional ABC transporters are composed of at least four subunits: two transmembrane (TM) domains where the transport process takes place and two cytoplasmic nucleotide binding domains (NBDs) providing the energy required for active transport (Figure 2). Although for a long time, all ABC transporters were assumed to share the same molecular mechanism, recently accumulating evidence suggests that the mechanisms might vary from one subfamily to another.[188−191] Despite the structural and likely mechanistic discrepancies among different ABC transporters, their NBDs are highly conserved and are likely responding to ATP binding and hydrolysis similarly. That is, as evidenced by many crystal structures,[192,193] the NBDs form a closed dimer upon ATP binding, which opens after ATP hydrolysis, a scheme commonly described in different mechanistic models.[192,194−196] The change in the dimerization state of the NBDs upon ATP binding and hydrolysis is the driving force for the transport as these events are coupled to different conformational states of the TM domains and their interconversion between the IF and OF states.

The conformational response of NBDs to ATP binding and hydrolysis has been best demonstrated by the crystal structures of the isolated NBDs from the maltose transporter (MalK).[197,198] The observed conformations for the NBDs within the context of a full transporter for several ABC transporters are also consistent with this picture.[4,88,199] MalK has been crystallized as an ATP-bound, closed dimer,[197] which appears to convert to an open dimer when Mg2+ is supplied to the crystallizing solution to allow for ATP hydrolysis.[198] The structure of the ADP-bound form of the MalK dimer is clearly in an open form. However, this structure does not establish whether dimer opening is an immediate result of ATP hydrolysis or requires the hydrolysis products, Pi, to be released from the NBDs. Furthermore, protein–nucleotide interactions during dimer opening cannot be derived from the static pictures provided by these crystal structures. Moreover, because there are two nucleotide binding sites in the NBDs, it has been debated whether ATP hydrolysis in both active sites is required for dimer opening and, if so, whether simultaneous or sequential hydrolysis events in the two binding sites constitute the mechanism. The dynamical view provided by atomistic MD simulations thus makes it an optimal tool for exploring the detail of the structural transitions constituting the posthydrolysis events.

To investigate the conformational changes after ATP hydrolysis, the bound ATP was hydrolyzed in silico to ADP and Pi at either or both of the binding sites in several MalK simulations.[200] This important biochemical reaction (ATP hydrolysis) can be conveniently modeled in MD simulations. After an initial equilibration phase of the ATP-bound system, we stopped the simulations and, using the positions of the three phosphate groups as a reference, converted ATP to ADP and Pi through breaking the bond between the γ- and β-phosphates. The proper number of hydrogens and oxygens were added to account for the hydrolytic nature of the reaction and to model the hydrolysis products in their relevant titration states. The simulations were then resumed with the new system. It was found that after the ATP is hydrolyzed, the nascent Pi moves quickly and significantly away from the β-phosphate of ADP, an event that can be attributed to the electrostatic repulsion between these two separate molecular species. Meanwhile, the ATP-associated Mg2+ cofactor moves between the produced ADP and Pi, so that it loses the coordination of a highly conserved glutamine residue (Q82) of the NBD. The nascent Pi, on the other hand, loses a critical hydrogen bond to a strictly conserved serine (S135) in the opposite MalK monomer (Figure 7A). Because the NBD dimer is held together largely by the interactions between the ATP and both of the two NBD monomers, the loss of a key H-bond at the dimer interface as a result of ATP hydrolysis destabilizes the dimer interface to a degree that eventually results in dimer opening, typically in 30–70 ns. The opening of the nucleotide binding sites at the dimeric interface can be clearly demonstrated by the distance between the two highly conserved ATP binding motifs, one from each NBD monomer (Figure 7B).

Figure 7

Conformational changes induced by ATP hydrolysis in the NBD of ABC transporters. (A) The local rearrangements of the nucleotide binding site induced by ATP hydrolysis, in particular disruption of a key hydrogen bond between the γ-phosphate and a serine residue (the thick dotted line in the left panel) at the dimer interface, eventually trigger the separation of the two monomers. The distance between the hydrogen bond donor and the acceptor in each binding site is shown at the right. (B) Global conformational change induced by ATP hydrolysis. In this particular simulation, the opening of the dimer is evident at both nucleotide binding sites. The distances between two nucleotide binding motifs at the dimer interface are recorded as an indicator of the degree of dimer opening.

As any other microscopic molecular event, dimer opening is a stochastic process; the opening time varies in different simulations, and any of the two nucleotide binding sites (not necessarily the site at which ATP has been hydrolyzed) may open first. Closer examination of the trajectories revealed the origin of the stochasticity: dimer opening requires simultaneous disruption of several H-bonds between the nucleotide and the opposite NBD.[200] As each H-bond exhibits an independent pattern of formation and breaking, which is facilitated by and also coupled to random collisions of the surrounding water molecules, simultaneous breaking of all H-bonds in each simulation is observed to happen on a different time scale.[48]

These simulations characterize the sequence of molecular events in the nucleotide binding sites following ATP hydrolysis that eventually lead to separation of the NBD monomers. In addition, they provide novel insight into the ambiguity about the stoichiometry between ATP hydrolysis and substrate transport: it is possible to induce NBD opening with only one ATP hydrolyzed in an NBD dimer. Because NBD opening results in the conformational changes in the TM domains, it is suggested that hydrolysis of only one of the two bound ATP molecules may be sufficient for interconversion of the transporter between the OF and IF states and thus to drive substrate translocation. This conclusion might explain the functionality of several ABC transporters with one degenerate nucleotide binding site,[201] in which the ATP binding capability of both sites is preserved but hydrolytic activity is present in only one of the sites.

The conversion of ATP to ADP and has also been used in the simulations of the isolated NBD dimers of other ABC transporters.[202,203] Despite the use of lengthy simulations on multiple replicates, NBD dimer opening was not observed in those studies,[202,203] possibly because of specific simulation conditions, including the choice of the protonation state of the nucleotide. On the other hand, NBD dimer opening has been captured with MD simulations when the effect of ATP hydrolysis was approximated by either complete removal of the nucleotide[204−206] or replacement of the nucleotide with only ADP (that is, without the inclusion of Pi).[208,207] A further separation between NBDs has also been demonstrated in simulations starting with a nucleotide-free, semiopen NBD dimer.[209]

The simulations described in this section exemplify some of the ATP-induced large-scale conformational changes in isolated NBDs of ABC transporters that have been successfully captured by equilibrium MD simulations. For an intact ABC transporter embedded in a membrane, it is much more difficult to capture similar conformational changes within the limited time scales of MD, because of the close coupling of the motions of the NBD and the TM domains, with the latter being further dampened by the membrane. As a result, NBD opening within the context of a full ABC transporter has been observed in only a small number of previous MD studies,[48,210,211] while the majority of such studies report no significant opening of the NBDs upon replacement of ATP with ADP.[38,47,99,212] Therefore, a full description of conformational changes underlying the complete transport cycle is still beyond the reach of conventional MD simulations, and one has to resort to more advanced simulation techniques (e.g., nonequilibrium methods) to investigate such structural transitions within currently feasible computational resources.

Concluding Remarks

The scope of molecular simulation as a biophysical technique has significantly expanded in recent years. At present, all areas of modern biomolecular sciences benefit from close integration of simulation technologies. In this review, using a number of membrane transporters, we have exemplified the application of these methodologies to the investigation of biomolecular dynamics and function. Molecular simulation techniques are highly useful when employed hand in hand with experiments. MD simulation techniques have been shown to be useful not only in explaining experimental observations but also in making predictions that can be tested experimentally. Furthermore, as in other research methods, the obtained results from simulation studies can be affected by the approximations and limitations of the method, and therefore, it is important to, when possible, verify and test the models and hypotheses generated by simulation studies through experiments. Given this, the examples in this review as well as many other studies demonstrate that because of the uniquely high spatial and temporal resolutions of MD methods, carefully designed simulations can provide insight that would be otherwise inaccessible. The detailed pictures of molecular events visualized by molecular simulations are by and large impossible to elucidate using currently available experimental techniques. As technological advancements continue at a rapid pace, simulation techniques are set to become a nearly essential analytical tool for biophysical investigations.