Nucleotide-Induced Conformational Dynamics in ABC Transporters from Structure-Based Coarse Grained Modeling

yet be determined, the model was then applied to predict nucleotide-induced conformational motions. Upon binding of ATP-mimicking ligands the structure changed from a conformation in which the nucleotide-binding domains formed an open shape, to a conformation in which they were found in tight contact, while, at the same time, a pronounced rotation of the transmembrane domains was observed. This ﬁnding is supported by normal mode analysis, and, comparison with structural data of the homologous vitamin B 12 transporter BtuCD suggests that the observed rotation mechanism may contribute a common functional aspect for this class of ABC transporters. Although in HmuuV noticeable rearrangement of essential transmembrane helices was detected, there are no indications from our simulations that ATP binding alone may facilitate propagation of substrate molecules in this transporter. Possible explanations are discussed in the light of currently debated transport scenarios of ABC transporters.


INTRODUCTION
ATP-binding cassette (ABC) transporters are molecular machines which are located in cell membranes and involved in the active transport of chemical substances [1,2], thus representing key players for the maintenance of cellular homeostatis. Owing to their enzymatic activity of binding ATP molecules and catalyzing their hydrolysis reaction, they are capable of converting chemical energy into internal mechanical motions which are used to implement the transmembrane exchange of substrates. ABC transporters can function either as importers which organize the efficient uptake of molecules, e.g., nutrients critical for intracellular processes [3,4], or as exporters, pumping lipids, metabolic products, toxins or drugs out of cells [5,6]. ABC transporters have also been evidenced to be clinically relevant since they can contribute to the resistance of cancer cells to drugs, to antibiotic tolerance, or cause diseases such as cystic fibrosis [7,8].
The common architecture of ABC transporters reveals two cytosolic nucleotide-binding domains (NBDs), the ABCs, and a pair of domains which are integrated into the membrane (the transmembrane domains, TMDs). Binding of intracellular ATP molecules and their hydrolysis inside the cassettes controls the conformation of the transmembrane domains which contain the substrate-specific translocation pathway, thus implementing the respective function of the transporter [9][10][11].
The structures of several ABC transporters have become available in the past and it has been shown that the NBDs are similar and constitute the engine core of the transporter, containing conserved sequence motifs which are involved in interactions with ATP molecules [12][13][14][15][16]. The TMDs, however, can considerably differ in their sequence and architecture despite all of them having a central pathway for the translocation of substrates. Moreover, since a few structures of transporters in the presence of nucleotides exist, valuable insights into possible conformational changes were gained and the basis for an understanding of their ATPdependent activity could be established [17][18][19][20][21]. In addition to static structural data, molecular dynamics (MD) simulations have been performed for the NBD components of various ABC transporters [22][23][24][25][26][27][28][29][30][31]. It is well established that binding of ATP causes closing of the NBDs and opening of them is obtained after ATP hydrolysis. Although MD simulations are capable of describing ATP-induced structural changes in the isolated NBD core of ABC transporters, the entire conformational dynamics underlying the transport cycle involves global motions on timescales which, due to their heavy demand on computational resources, are still beyond the reach of such methods [32,33].
Therefore, despite extensive experimental and modeling studies, important functional aspects of the operation of ABC transporters are still lacking to date. This particularly refers to the coupling of ATP-induced activity of the NBDs to the transport activity inside the TMDs. Resolving the conformational motions underlying these processes is essential to understand the mechanism of ABC transporters. In this situation, modeling studies that are based on present structures of transporter proteins and provide approximate descriptions of their dynamics may play an important role.
To overcome the apparent difficulties associated with MD methods, coarse-grained descriptions of protein dynamics have become more and more popular within the last decades [34,35]. A widespread approach to analyze conformational motions in proteins with feasible computational burden is the elastic network model (ENM) [36,37]. In this approximate description, entire amino acids are typically modeled as single point-like particles and the complex intramolecular interactions are replaced by empiric harmonical potentials which mediate effective interactions between them. Thus, within the framework of this mechanical model, a protein is viewed as a network of beads connected via deformable elastic springs. It has been shown that ENMs are able to reproduce well the pattern of residue displacement in protein structures due to thermal fluctuations and describe ligand-induced collective dynamics in molecular machines [38][39][40]. In most former studies predictions of functional conformational motions in proteins were based on the computation of normal modes of elastic networks, an approximation which allows for an efficient numerical implementation, but nonetheless may also be prone to fail as reported earlier [41,42]. Elastic network normal mode analysis has been previously applied to study nucleotide-dependent conformational dynamics of the vitamin B 12 ABC transporter BtuCD [43]. More recently, relaxational elastic network models which take into account the full non-linear network dynamics have been applied to study mechanochemical motions in protein machines [44][45][46]. In a recent study we have used such models to follow, for the first time, entire operation cycles of an important motor protein in a structurally resolved manner [47,48].
Here we propose and apply an approach which, based on the relaxational elastic network model, intends to investigate nucleotide-dependent dynamics of ABC transporters at the coarse-grained level. Based on the crystal structure of the nucleotide-free state of the transporter we have constructed a ligand-elastic-network model that implements interactions with ATP molecules in a simplified fashion. The focus in the performed computer experiments was on the study of conformational motions induced by binding of ATP ligands and the subsequent response generated in the transmembrane domains. By tracing the formation of an approximate in silico version of the nucleotide-complexed state of an ABC transporter the proposed approach allows to follow conformational changes of the full protein structure and thus explore the coupling of nucleotide-induced dynamics in the NBDs to motions inside the TMDs and relate it to possible gating mechanisms. As a proof of concept the model was first applied to transporter structures for which multiple structures are known and the quality of its predictions was evaluated by comparing the ligandcomplexed conformations obtained from the simulations with the crystal structures in the presence of ATP-analogs. In the next step we have performed simulations for the ABC protein HmuUV, a transporter involved in conducting the import of heme substrates, which is essential in bacteria [49], but for which yet only the nucleotide-free structure could experimentally be determined [50].
Particularly our findings summarize as follows: (i) generally, binding of nucleotide-mimicking ligands induced largeamplitude global and ordered conformational motions that were robust, i.e., they were not sensitive to the specific interaction pattern applied between ligands and the protein within our description, (ii) model predictions agree remarkably well with structural data; for the MalK ABC dimer the switching from open to closed conformation is reproduced and for the fulllength maltose transporter the progression from inward-facing to outward facing state of the transmembrane domains is correctly obtained, (iii) in HmuUV a coupling of tweezer-like motion of the NBDs and rotational motions of the TMDs was observed upon binding of ligands; this mechanism of coupled domain dynamics is supported by normal mode analysis and comparison with available structures of the homologous vitamin B 12 transporter BtuCD, (iv) in HmuuV noticeable rearrangement of essential transmembrane helices was detected but there are no clear indications from our modeling of how conformational motions induced by ATP binding alone may relate to mechanical aspects of gating in this transporter.

Network Construction and Dynamics
The performed investigations of ABC transporters were based on the crystal structure of the respective nucleotide-free conformation. We describe here the construction of their corresponding elastic networks and refer to Table 1 for more detailed information. The coarse-grained elastic network of a protein was obtained by first replacing each amino acid residue in the crystal structure by a single bead which was placed at the position of the alpha-carbon atom of the respective residue. These equilibrium positions were denoted by R (0) i for bead i. Then, to determine the pattern of network connections, the Network beads 734 1351 1108 Summary of the studied ABC transporter structures with their corresponding PDB ID code and the residues considered in the network construction. Residues from conserved motifs used in the modellng are provided.
j | between equilibrium positions of any two beads i and j was compared with a prescribed interaction radius r int . If this distance was below r int , the two beads were connected by a deformable elastic spring. The network connectivity was stored in matrix A with entries A ij = 1, if beads i and j are connected, and A ij = 0, else. For all ABC protein networks we have used an interaction radius of 9Å. The choice of the interaction radius in elastic network simulations has been widely debated on in the field. If the value is too small, the network would decompose into lose or even disconnected substructures, which would lead to unreliable results. On the other side, a too large interaction radius would render the network very stiff, thus preventing functional motions to occur. In this work, the interaction radius has be chosen based on observations from a previous study, where values between 8Å and 9.5Å have been found to be optimal to describe domain motions in many proteins [51].
The total elastic energy of the network is the sum of all spring contributions, i.e., U = Here, N is the number of beads in the protein network, d ij = | R i − R j | is the actual length of a spring connecting beads i and j in some deformed network conformation, with R i being the actual position vector of bead i, and d (0) ij is the corresponding natural length (as extracted from the PDB file). The spring stiffness constant κ is assumed to be the same for all network springs.
When thermal fluctuations and hydrodynamical interactions are neglected, the dynamics of the network can be described by a set of over-damped Newton equations. For bead i the equation of motion is On the left hand side of Equation (1), γ is the friction coefficient assumed to be equal for all network beads. The first part on the right hand side are the elastic forces exerted by network springs which are connected to bead i and the second term F i,lig is the force exerted by an ATP ligand. The parameter σ i equals 1 if a ligand is connected to bead i, and is 0 otherwise. Interactions between the protein elastic network and ligands are described in the following section. Explicitly, Equation (1) reads

Modeling of Nucleotide Binding
Based on the constructed elastic network of a nucleotide-free state of the ABC transporter we have established a ligand-elasticnetwork model in which each of the two ATP molecules was described by a single ligand bead and its binding process was modeled by placing it into the center of the nucleotide binding pocket and connecting it to pocket beads via elastic springs. The binding configuration is schematically displayed in Figure 1. In the simulations, ligand bead 1 was connected to the Walker A motif beads (Gly-Pro-X-Gly-X-Gly-Lys-Ser) from one NBD and the ABC motif beads (Leu-Ser-Gly-Gly-X-X-Gln-Arg) from the other domain. Accordingly, ligand bead 2 was connected to the corresponding motif beads from the second pocket. The natural lengths of the additional ligand-springs were chosen smaller than their corresponding initial lengths, implying that initially they were stretched and upon binding the generated forces between each ligand and the selected pocket beads were attractive.
To check the robustness of our ligand-model, we have followed for each investigated ABC transporter, the formation of the ligand-complexed structure for different initial ligand-binding conditions. Starting always from the native conformation of the ABC protein network, the two ligand beads were simultaneously placed in their pockets but the natural lengths of ligand springs have been randomly varied in each realization (i.e., the strengths of initial deformations of ligand springs have been altered). More precisely, the following scheme was implemented in the simulations. The natural lengths of the ligand springs have been chosen as d (0) to the motif bead with index p 1 k (p 2 k ) upon its placement. The parameters α p 1 k and α p 2 k were random numbers chosen from an interval between 2.0 and 5.0, ensuring that natural spring lengths were shorter than their initial lengths. To impose symmetric ligand interactions, these parameters were the same for both ligands (i.e., α p 1 k = α p 2 k ), meaning that the spring connecting ligand 1 to a particular motif bead and the corresponding spring connection for ligand 2 had roughly the same prescribed natural length. Then, as a normalization condition for all realizations, the natural spring lengths were rescaled, such that the amplitude of the total force F ini lig exerted by each ligand upon binding takes a fixed common value, i.e., | F ini The equation for ligand 2 was the same with index 2 instead of 1. Here, R lig 1 is the actual position of ligand 1 and d lig 1 ,p 1 | are the actual lengths of ligand springs. For the ligand beads the same friction coefficient γ as for other network beads was assumed. The additional springs also had the same stiffness constant κ as all other elastic springs. According to Newton's actio and reactio principle the forces exerted by the ligand in Equation (2) are related to the forces in Equation (3) (2) and (3) the dependencies on the friction coefficient γ and the spring stiffness constant κ can be removed using an appropriate rescaling of time.
Subsequent to binding of the two ligand beads, we have followed the formation of the in silico ligand-bound complex of the ABC transporter by integrating the set of Equations (2) and (3) until a stationary structure of the network-ligand complex was reached (time moment 2500 in the simulations). This procedure was carried out 50 times (for MalK and MalFGK 2 proof of concept simulations 20 times, respectively), each time with a different set of coefficients α p 1 k . In the simulations the equations of motion were numerically solved by employing an first-order integration scheme with a time-step of 0.1.
To visualize conformational changes resulting from ligandbinding, one ligand-bound complex obtained from a single simulation in each investigated ABC transporter was selected and its structure shown in the ribbon representation. To ensure reproducibility of these results, the parameters of ligand interactions used in the specific simulations are provided in Supplementary Table 1.

Characterization of Conformational Motions
In the simulations conformational motions have been characterized by tracing distance changes between network beads (or groups of beads). To determine the size of a nucleotide binding pocket, its gyration radius was computed as the RMSD of network beads which corresponded to the Walker A and ABC signature residues of a pocket (see Table 1) with respect to its geometrical center. The separation between the two NBDs was measured as the distance between their geometric centers. In the MalFGK 2 maltose transporter the width of the cytosolic gate was determined as the distance between the two network beads which corresponded to residues Thr176 from chain G and Leu429 from chain F, respectively. For HmuUV several variables have been monitored during the simulations. Beside the NBD separation we have also measured the separation between the TMD coupling helices as the distance between the Asp222 network bead of each TMD. The rotational shift of the TMDs with respect to the NBDs has been quantified by the angle between the vector which connects the geometric centers of the NBDs and that which connects the two centers of the TMDs. To further describe global motions of the TMDs, we have calculated during the simulations the principal axes of each TMD as the eigenvectors of its inertia matrix. To detect motions localized at the central interface of the TMDs, the distance between the Ala146 residue beads which belong to the TM5 helices and are located in the vicinity of the putative cytoplasmic substrate gate has been monitored. Near the periplasm side the distance between the Thr167 residue beads (from the TM5 helices) and that between the Arg176 beads (from the H5a helices) was traced. To depict conformations obtained from the simulations in the ribbon representation, a structural reconstruction from the alpha-carbon trace was performed [52]. For the visualization the VMD software was used [53].

Normal Mode Analysis and Modified ENM
The normal mode analysis for the HmuUV structure was performed using the AD-ENM server. There, the same interaction radius as in our simulations was used for the elastic network and the amplitude of displacement was set to 3.0. The structural displacement of the lowest non-zero mode was compared with the findings obtained from our simulations.
For the HmuUV ABC transporter, in addition to the usual elastic network model, a modified setup was considered. In the conventional elastic network description the bead connectivity is determined by applying the interaction radius for a specific static conformation of the protein and, once fixed, it is assumed not to change under structural changes of the network, i.e., existing links between beads cannot be broken and new links cannot be established either. Hence, the separation of two network beads which are initially connected by a spring is bound by the harmonic potential coupling them and can never be too large. In the elastic network of the HmuUV apo structure several links connecting the two TM5 helices at the central TMD interface close to the putative cytoplasmic gate were present. Since such links could artificially stiffen the network in that region and, hence, possibly prevent relative motions of the transmembrane helices, a modified elastic network description in which the interface between the two TMDs was softened was considered. In this model, for simplicity, the few links which connect the TM5 helix from one TMD with that of the other TMD have been cut after the original construction of the elastic network. More precisely, the 15 links which initially connected residue beads Asn144 to Ile153 of one TMD to the corresponding residue beads in the other TMD were removed. Then, in an independent set of simulations, the analysis for HmuUV was repeated using the modified description. Figure 1 shows the common architecture of ABC transporters consisting of two nucleotide-binding domains (NBDs) located in the cytosol and a pair of transmembrane domains (TMDs). In each NBD the conserved Walker A and ABC signature residue motifs, which are critical for the binding of ATP molecules and their hydrolysis reaction can be found. In the absence of nucleotides the motifs are spatially separated and the two binding pockets formed by them are in an open conformation, each ready to bind an ATP molecule. The TMDs typically consist of several alpha helices. At the central domain interface some of them constitute a confined channel structure through which the substrate is transported across the membrane.

RESULTS
We have considered the elastic network description to investigate conformational dynamics of ABC transporters in computer experiments. The network consisted of identical beads which represented entire amino acid residues and elastic springs connecting any two beads if their spatial distance was below a prescribed interaction radius. Details of the network construction can be found in the Materials and Methods section. The elastic network model is an approximative mechanical description of protein dynamics, capable of computationally handling large proteins with several hundreds (or even thousands) of amino acids and following slow conformational motions inside them in a structurally resolved way. The internal motions of proteins in solution are characterized by low Reynolds numbers and therefore dominated by friction. In fact, it has been shown that inertial effects are negligible on time scales of tens of picoseconds and slower conformational motions become relaxational [54]. Therefore, the slow conformational motions with characteristic timescales of milliseconds taking place in protein machines should be purely damped. Neglecting hydrodynamical effects and thermal fluctuations, the dynamics of the protein network can therefore be described by mechanical Newton's equations of motion considered in the over-damped limit. Thus, within the approximate nature of the model, internal motions in proteins are described as processes of conformational relaxation of its corresponding elastic network.
The aim of this study was to investigate nucleotide-induced conformational dynamics of ABC transporters, analyse the coupling of domain motions and discuss how structural changes inside the transmembrane part may contribute to the mechanism of substrate transport. For that purpose, we have constructed a ligand-elastic-network model of ABC transporters which incorporated binding of two ATP mimicking ligands to the pockets inside the NBDs. This approach will be described in the following section.

ATP-ligand Binding and Related Conformational Motions
Due to, in part, the lack of experimental data, interactions with ATP and the effect of binding and its hydrolysis is hitherto unknown in many ABC transporters. However, the two nucleotide-binding pockets are well-defined by the presence of the conserved Walker A and ABC signature motifs located in each NBD. Moreover, it has been shown previously, that most ABC transporters bind two ATP molecules in each single catalytic cycle [55]. On the other side, the coarse-grained nature of the elastic network model does not allow to resolve the complicated interactions between the protein and ATP molecules in full detail. To investigate conformational motions induced by binding of ATP molecules, we have constructed a dynamical ligand-elasticnetwork model of ABC transporters which incorporated binding of two ATP nucleotides to the NBDs in a simplified fashion. Details of the modeling and the network dynamics are provided in the Materials and Methods section and the limitations of our approach are mentioned in the Discussion section.
Similar to our previous publication, where entire operation cycles of the hepatitis C virus helicase motor protein were investigated using ENM modeling [47], binding of ATP was phenomenologically accounted for in the present study. Consistent with the coarse-grained character of the elastic network model, we have described an ATP-ligand by a single bead and modeled its binding process to the nucleotide pocket by allowing it to physically interact with network beads present there. We have selected those beads as interaction spots that corresponded to residues from the conserved sequence motifs mentioned above. More precisely, a ligand could interact with beads of the Walker A motif from one NBD and beads of the ABC motif from the other NBD. Binding of the two ATP-ligands was implemented by simultaneously placing each ligand bead into its binding pocket and creating spring connections to the motif beads. These additional elastic springs were initially stretched so that upon binding, attractive forces between each ligand and the selected pocket beads were generated. Therefore, as a result of the binding process, network deformations were initiated. First localized in the two binding pockets, they gradually spread and conformational relaxation motions of the entire protein network took place until a stationary structure of the elastic network complexed with the two ligands was reached.
To follow the formation of an approximate in silico version of the ligand-complexed transporter structure upon binding of the two ligand beads was the focus of our computer experiments.
We have numerically integrated the set of equations of motion to obtain the position of all network beads at every instant of time. Since in our approach interactions between ATP ligands and the pocket were modeled only roughly, we have tested the robustness of the ligand network description and its predicted conformational changes by following the formation of the ligand-bound complex many times, each time starting from the nucleotide-free network. In each realization different ligand-binding patterns were applied by randomly varying initial deformations of the ligand springs (see Materials and Methods). To characterize conformational motions subsequent to binding of the two ligand-beads we have traced several geometric variables during the simulations, which reflect both global domain rearrangements as well as structural changes localized in the substrate translocation cavity of the TMDs (see Materials and Methods). A single process of conformational relaxation in the transporter elastic network, which was induced by binding of the two ligand beads and corresponded to the formation of a transporter-ligand complex, could therefore be visualized by a single trajectory in the space of these variables.

Proof of Concept Simulations
Predictions of the proposed ligand-elastic-network model were benchmarked by applying it to ABC transporters for which several crystal structures are known. Conformational motions subsequent to ligand binding were analyzed and the structure of the ligand-bound complex obtained from the simulations was compared with that of the corresponding crystal structure determined with ATP-analogs. First, for a general test of the suggested method, the isolated ATPase subunit of an ABC transporter was considered. After that, the full-length maltose transporter MalFGK 2 was investigated. A summary of the studied structures and details of the simulations are given in Table 1.

MalK ATPase Subunit
The MalK homodimer represents the ATPase engine which powers substrate translocation in the maltose transporter system of Escherichia coli. Each MalK monomer consists of a nucleotide binding domain which harbors the Walker A and ABC signature motifs and a regulatory domain (RG) which interacts with other enzymes [13]. Crystal structures of MalK in the absence of ATP and complexed with nucleotides are shown in Figure 2. Without ATP the monomers are separated and the MalK dimer structure is in the open shape conformation (Figures 2A,B). When binding of two ligands in the elastic network which corresponded to that open structure was modeled, we observed for all different realizations of ligand-binding conditions similar conformational dynamics along the formation of the ligandcomplexed structures, as evidenced by the low dispersion in the trajectories which monitor changes of the protein structure (shown in Figure 3). Subsequent to binding of the ligands, deformations of the elastic networks which were localized in the nucleotide binding pockets became generated and the pockets closed. We detected a decrease in the gyration radii of the pockets of about 6Å (see Figure 3A). When such initial deformations eventually spread over the entire MalK network they induced global large-amplitude relative motions of the two NBDs with the separation between them decreasing about 9Å (see Figure 3B). Therefore, in all final ligand-complexed networks the MalK dimer was in a closed conformation with the two NBDs found to be in tight contact. All predicted ligand-bound structures obtained from the simulations compared with the ATP-bound form of the MalK crystal structure by a RMSD value less than 4Å; in comparison, the initial ATP-free MalK and ATP-bound MalK differ by an RMSD of 7.2Å. In Figures 2a,b we show one selected in silico ligand-bound structure together with the ATP-complexed crystal structure. The two conformations agree remarkably well. They align with an RMSD of 2.7Å. In the in silico structure the nucleotide binding pocket have gyration radii of 6.4Å and 6.3Å, respectively, and the NBDs are separated by 20.4Å. For comparison, the corresponding values in the ATPbound crystal structure are 6.7Å and 6.8Å for the gyration radii of the pockets and 22.4Å for the NBD distance.

Maltose ABC Transporter
In the next step of benchmark simulations we have checked how well the coupling of ATP-induced motions and transmembrane dynamics is reproduced in our model by applying it to a fulllength ABC protein. We have chosen the transporter MalFGK 2 from Escherichia coli which is involved in the import of maltose/maltodextrin substrates and represents the probably best studied ABC systems, for which multiple structures are available [20]. Its overall structure comprises the two ATPase forming cytoplasmic MalK subunits and two transmembrane domains, MalF and MalG. The structures we used have been co-determined with the periplasmic maltose binding protein MalE. Since here we are not interested in a detailed explanation of the transport cycle we have omitted the periplasmic parts in the simulation. Maltose importer structures are shown in Figure 4. Without nucleotides the MalK dimer is in the expected open conformation (Figure 4A). At the same time the TMDs are arranged in an inward-facing conformation with a closed periplasmic gate and the gate to the cytoplasm being wide open (Figure 4B).
When the nucleotide binding process was modeled in this transporter, we found similar conformational dynamics upon binding of the two ligands as for the previous simulations of the isolated MalK dimer; the nucleotide binding pockets were shrinking in its size and the large scale global motions which bring the NBDs in tight contact set in. In Figures 5A,B the radii of the pockets and the distance between the NBDs are plotted as they change along the formation of the ligand-bound complex. Again, as in the previous simulation, the conformational changes were robust and not sensitive to variations in the ligand binding interactions in our model. The gyration radii of the pockets changed about 2.5Å and in ligand-complexed conformations the pocket sizes compared well with those found in the ATPcomplexed crystal structure (see Figure 5A). The separation of the NBDs decreased by approximately 2.5Å to a value around 22Å in the ligand-complexed conformations (see Figure 5B), which is in well agreement with the NBD distance computed for the ATP-complexed crystal structure (22.3Å). Regarding the coupling of nucleotide-induced motions in the NBDs to motions inside the TMDs, we found that all predicted ligandbound structures obtained from the simulations compared with the ATP-bound form of the maltose transporter by a RMSD value of less than 3Å, indicating that the model reproduced well the nucleotide-activated conformational rearrangements of the full transporter structure. To check details of the transmembrane gating we have monitored in our simulations the size of the cytoplasmic gate along the formation of the ligandbound transporter structures and found that it decreased by approximately 1.5Å (see Figure 5C) and is centred around a value of 8Å in the ligand-complexed states. For comparison, the gate size in the ATP-complexed crystal structure is 7.6Å.
The obtained results evidence that our model is able to reproduce the ATP-induced structural changes in the maltose transporter system remarkably well. To visualize the structural agreement we have selected a single in silico ligand-bound complex of the transporter and compared its conformation with that of the ATP-complex crystal structure (see Figure 4a). Overall the conformations align with an RMSD value of 2.6Å. Figure 4a illustrates the high degree of conformational similarity of the transporter structures. This applies to the NBDs where the ligand binding events take place as well as to the remote TMDs, where the model does correctly reproduce the progression from the inward-facing to the outward-facing channel shape (compare Figure 4B and Figure 4b). As for the quantitative agreement we found that in the in silico structure the nucleotide binding pocket have gyration radii of 6.8Å and 6.7Å, respectively, and the NBDs are separated by 22.0Å. For comparison, the corresponding values in the ATP-bound crystal structure are 6.7Å and 6.8Å for the gyration radii of the pockets and 22.3Å for the NBD distance. The size of the cytoplasmic gate was 7.9Å (in silico) and 7.6Å (crystal).

HmuUV Heme Transporter
The performed benchmark simulations have evidenced the quality of the model to predict functional conformational dynamics of ABC transporters remarkably well. Next, the model is applied to the HmuUV transporter for which only a single structure, that of the nucleotide-free state, has yet been determined [50]. HmuUV is an ABC importer which is expressed in bacteria as part of their survival strategy since it presents an efficient agent for the acquisition of iron, which typically has a low intracellular concentration but is needed for various chemical processes [49]. Figure 6A shows the domain architecture of the HmuUV heme transporter in its nucleotide-free conformation. In the NBDs (chains C and D) the conserved ATP binding motifs (Walker A and ABC signature motif) are spatially separated and form two open binding pockets for ATP ligands at the domain interface (see Figure 6C). At the interface of the two TMDs (chains A and B), transmembrane helices 5 (TM5) and helices 5a (H5a) form a cone-shaped cavity which has been assigned the putative heme translocation pathway [50]. Their outwardfacing arrangement indicates that in the nucleotide-free state of the transporter the pathway is open to the periplasm, whereas it is sealed to the cytoplasm by a narrow gate and heme passage would thus be sterically hindered in this region ( Figure 6B).
Like in the previous simulations, following the construction of the HmuUV elastic network, ligand binding to the NBDs was modeled and the subsequent conformational changes followed by monitoring relaxation trajectories. The results are shown in Figure 7 and again, like in the previous cases, they were not sensitive to variations in the ligand binding conditions.

Coupled Domain Motions
As for the motions generated inside the NBDs, they resemble those found for the MalK and the maltose transporter. The size of both nucleotide binding pockets decreased about 4Å (see Figure 7B) and large-amplitude relative motions of the NBDs bringing them into tight contact, with the distance between them decreasing about 8Å, were observed (see Figure 7A). Considering the coupling of domain dynamics in HmuUV, we find that ligand-induced structural changes in the NBDs cause drastic consequences in the conformation of the remote TMDs. First, the distance between two coupling helices which each connect a TMD with its NBD was significantly decreased about 6Å (see Figure 7C). Second, at the same time as the NBDs approached one another, the angle between the NBDs and the TMDs changed by about 20 • , indicating global rotational motions of the TMDs relative to the NBDs (see Figure 7A). Rotational motions are first small after ligand binding but become pronounced when the steady conformations of the ligand-bound complexes are approached. In addition to the rotational dynamics, tilting motions of each individual TMD were observed as evidenced by a change in the orientation of the long principal axes of ∼ 5 • (see Figure 7D). To visualize global conformational changes, the structure of one selected ligandbound complex was superimposed with that of the apo crystal structure of the HmuUV transporter (shown in Figure 8A). In the in silico obtained ligand-bound conformation the NBDs are separated by 27.3Å and the coupling helices are 31.3Å apart. The corresponding distances in the apo structure are 35.3Å and 37.8Å, respectively. Further on, the structural comparison clearly evidences the coupling of linear tweezer-like motions of the NBDs and a rotational displacement of the TMDs (see Figure 8a). This mode of coupling is confirmed by the normal mode analysis (see Materials and Methods). Indeed, the conformational motions corresponding to the first normal mode of the HmuUV elastic network (visualized in Supporting Video S1) resemble the collective domain motions identified by our model. Moreover, for the vitamin B 12 ABC transporter BtuCD, which is structurally similar to HmuUV, a similar scheme of domain coupling can be observed, as revealed through comparison of apo and nucleotide-complexed conformations [56,57] (presented in Figures 8B,b). Such aspects will be further commented on in the Discussion section.

Rearrangement of TM Helices
In order to detect possibly functional motions localized in the putative heme translocation cavity, distances between residues near the periplasmic and cytoplasmic gate were traced as they changed along the formation of the ligand-complexed HmuUV structures. No significant change could be detected in the distance between the Ala146 residues of the TMDs, which are in proximity of the cytoplasmic gate (see Supplementary Figure 1A). Similarly, variations in the distance between residues Arg176 and Thr167, both located near the exit to the periplasm, were small (about 1Å see Supplementary Figure 1B). In all ligand-complexed structures latter distances were decreased as compared to the apo structure, indicating that relative motions of the TM5 and H5a helices may tend to cause slight closing near the periplasm side. The structural rearrangement of TM helices is shown in Figure 9, where the TMDs of the HmuUV apo structure was superimposed with those of the single ligandbound structure obtained from the simulation (the same ligandcomplexed structure as shown in Figure 8). Coupled to ligandinduced conformational dynamics of the NBDs, the internal organization of all TM helices and their relative orientations was changed. As shown in Figure 9, for TM helices 1-4 and 6-9 the changes were pronounced, whereas at the central domain interface helices 5, 5a, and 10 underwent only marginal tilting and no significant rearrangements inside the putative heme translocation cavity could be detected. These observations are confirmed by the normal mode analysis (see Materials and Methods). The structural displacements inside the TMDs, as predicted by the lowest normal mode, is found to agree remarkably well with results from our ligand-elastic-network simulations (compare Figure 9 left and middle panel).
In the original elastic network description, which the performed ligand-binding model was based on, the two TM5 helices were connected by several elastic springs in the region close to the putative cytoplasmic heme gate. Since such links could possibly prevent functional relative motions of TM helices, a modified elastic network description in which the interface between the TMDs was softened was considered. Details of the modified model are given in the Materials and Methods section. In an independent set of simulations, the formation of ligandcomplexed HmuUV structures was followed using the modified description. Trajectories displaying changes in the widths of the cytoplasmic and periplasmic gate during the new simulations are provided in Supplementary Figures 1C,D. Despite the softened TMD interface they are still small (around 1Å) and their significance in functional gating cannot be evidenced. Structural changes inside the TMDs, which correspond to a single ligandbound conformation obtained from the modified ENM, are shown in Figure 9 (right panel).
In conclusion, rearrangements of TM helices as obtained from the undertaken modeling do not give clear evidence for either gate opening or closing at the cytoplasm and periplasm side and hence, results from the current analysis may not be able to provide an explanation of how nucleotide-dependent conformational motions are coupled to the transmembrane gating dynamics in the HmuUV ABC transporter (see Discussion).

DISCUSSION
In this paper a coarse grained model which allows to investigate nucleotide-induced conformational dynamics in ABC transporters is presented. Based on the elastic network description of the apo structure of the transporter, binding of nucleotide ligands to the ATP binding domains was mimicked in an approximate fashion. In this way the formation of the ligand-bound transporter complex could be followed and the underlying conformational motions analysed in computer simulations. The explanatory quality of the proposed model was first confirmed in a set of simulations where nucleotideinduced conformational changes predicted by the experimentally obtained crystal structures of the MalK NBD dimer and the full-length maltose ABC transporter were reproduced extremely well in silico. The model was then applied to explore functional dynamics of the heme-specific ABC transporter HmuUV, for which yet only the apo structure could been determined. It is remarkable that despite the gross simplifications present in the elastic network model, i.e., a protein is viewed as a  meshwork of beads connected by deformable springs, it can still well describe the ATP-induced functional conformational changes in ABC transporters. For the MalK subunit the switching from open to the closed dimer configuration upon ATP binding is reproduced. Moreover, for the full-length maltose transporter MalFGK 2 progression from inward-facing to the outward-facing state is correctly obtained. Hence, the model is apparently also capable of covering functionally important aspects of intramolecular domain communication and coupling of nucleotide-induced conformational motions in the NBDs to gating dynamics in the remote TMDs in ABC transporters. As a common aspect of the conformational dynamics we found that in all ABC simulations binding of nucleotide ligands induced large-amplitude conformational changes inside the transporter structure which corresponded to ordered domain motions. These motions were robust and not sensitive to the specific interaction pattern that was applied between ligands and the transporter within our description. This indicates that interactions with nucleotides generate a generic response, a property which may assure that similar well-organized domain motions are repeated in each ATP-induced operation cycle and substrate transport is robustly maintained, despite the presence of possible perturbations. The most significant response was identified inside the two NBDs which, subsequent to ligand binding, underwent substantial motions bringing them from an open-shape conformation in which they were well separated, to a closed-shape form where they were found in tight contact. This structural transition is essential for the function; only in the compacted dimer formation with the nucleotides being tightly sandwiched between the Walker A and the ABC signature motifs, their hydrolysis reaction can be performed. Closing/opening of the NBD dimer related to ATP binding/hydrolysis represents a well-established mechanism underlying the cyclic operation of ABC transporters. However, the coupling to conformational dynamics in the transmembrane domains and their role in substrate gating is generally explored far less.
For the maltose transporter MalFGK 2 which represents the probably best studied ABC system with multiple structures being available, our model correctly reproduced major aspects of the conformational coupling: in the absence of ATP ligands the NBDs are in the open state and the TMDs are in the inward-facing conformation where the substrate may leave the transporter through the opened cytoplasmic gate, and, at the same time, the gate to the periplasm is closed and substrate entrance from the extracellular space is permitted. In contrast, in the nucleotide-complexed transporter the NBDs are found  Figure 8, obtained from the usual ENM. Structure changes predicted by the first normal mode are shown in the middle column. Right column: A single ligand-bound conformation as obtained from the modified ENM. Three different perspectives are provided. In row (A), the front view is chosen and the distance between the coupling helices is given. The view from the periplasm and that from the cytoplasmic side are provided in rows (B,C), respectively (for clarity the coupling helices have been omitted in these representation). For one TMD the TM helices are assigned their numeric labels.
in close proximity, but the TMDs underwent a transition upon which the gate to the cytoplasm became closed and essential transmembrane helices were tilted toward an outward-facing conformation, possibly opening the gate to the periplasm.
For the heme ABC transporter HmuUV the situation is quite different since experimental data is basically lacking and only the crystal structure in absence of ATP analogs is known yet. In this structure the NBD dimer is opened and the TMDs adopt an outward facing conformation where the putative substrate translocation cavity is accessible from the periplasm but closed to the cytoplasm. Our simulations reveal an interesting coupling mechanism of domain dynamics in which the linear tweezerlike closing motions of the two NBDs, induced by nucleotide binding, are translated into pronounced global rotation motions of the TMDs. By comparing available structural data we deduced a similar coupling scheme for the vitamin B 12 transporter BtuCD, indicating a possible communality of domain communication in this class of ABC transporters. Rotation of the TMDs with respect to the NBDs should generally generate twisting at their interface and possibly alter the orientation of transmembrane helices. In the HmuUV simulations tilting of the TMDs was found, the coupling helices underwent significant motions, and structural rearrangements of TM helices were observed. Nonetheless, the conformation of the putative translocation cavity was not much affected by these motions. Hence, the current analysis is not able to provide indications that may link nucleotide binding to substrate gating activity in HmuUV.

Transport Mechanism in ABC Importers
Although a complete model of ABC transporters has not yet been fully elucidated, the transmembrane exchange of substrates can be described on a general level in terms of the alternating access model. It assumes the conformation of the substrate translocation cavity to exist in two states: In the outward-facing state, substrate molecules may enter the TMDs through the accessible periplasmic side, whereas the gate to the cytoplasm is sealed. In the inward-facing state instead, the gate to the perisplasm is closed but the cavity to the cytoplasm is opened, allowing the substrate to leave into the intracellular space. Substrate translocation is controlled by the switching between both states, a process which is driven by the nucleotide-dependent conformational dynamics of the NBDs, and hence couples to the open/closed state of the NBD dimer. For the transport cycle of ATP importers two opposing coupling mechanisms have been established based on structural and biochemical data for the maltose and the vitamin B 12 transporter. Type I importers (e.g., the maltose transporter MalFGK 2 ) adopt the inwardfacing conformation in the absence of nucleotides but are found switched to the outward-facing conformation in the ATP-bound state [13,58]. A distinct mechanism is debated for type II importers such as BtuCD [21,56]. In BtuCD the outward-facing conformation is revealed in the apo state, in which the substrate cavity at the cytoplasmic side is sealed by gate I, but a second gate (gate II) is open [57]. In the AMP-PNP-bound state, the outward-facing conformation is obtained likewise but this time the space to the cytoplasm is closed by gate II whereas gate I became opened. Therefore, in the framework of the BtuCD model the process of nucleotide binding is neither coupled to gating activity near the cytoplasm, nor does it correlate with significant changes near the periplasmic region. As a matter of fact, for BtuCD there is currently no data available which evidences the release step of the vitamin B 12 substrate through the cytoplasmic gate of the transporter.
Similar to BtuCD, in the heme transporter HmuUV the outward-facing conformation is found coupled to the apo opendimer state of the NBDs. With respect to the substrate gating activity our results for HmuUV are consistent with the BtuCD model, i.e., nucleotide-binding does not significantly affect the conformation of TMDs near the periplasm and changes near the cytoplasm cannot be directly related to substrate gating. It should be noted that for the conformational dynamics in HmuUV, the role of the substrate binding protein HmuT, which delivers the heme molecules from the extracellular space by docking at the periplasmic site of the TMDs, is unknown at present [50]. In the transport cycle of the BtuCD model the corresponding substrate binding protein BtuF is assumed to trigger rearrangements of the TMD translocation cavity [56,57].
Results obtained here for the HmuUV ABC transporter have predictive character and experimental verification is required. Certainly it would be important to determine the structure complexed with ATP-analogs and in the presence of HmuT. Furthermore, predictions of domain interplay, in particular the rotational dynamics, can possibly be checked in site-directed spin labeling electron paramagnetic resonance (EPR) experiments.

Model Limitations and Concluding Remarks
Any coarse-grained model is obtained through approximations and therefore underlies limitations. The elastic-network model provides a mechanical description of nucleotide-induced dynamics in protein machines whereas, as the price paid for of the applied simplifications, the involved chemical processes cannot be resolved. In this study the elastic-network model could not resolve the chemical details that are important to correctly model interactions of ABC transporters with ATP molecules. Therefore, similar to our previous elastic-network study where full operation cycles of the hepatitis C virus helicase protein motor could be modeled [47], the process of ATP-binding was emulated only roughly and physical interactions were accounted for by empiric effective potentials. On the other hand the obtained results suggest that details of the ligand-protein interactions are not particularly important when considering nucleotide-induced global conformational dynamics; the same robust and generic motions were generated and they were not sensitive to variations in the interaction pattern. This observation agrees with our elastic-network studies of mechano-chemical dynamics of other molecular machines [45,47]. It is indeed remarkable that a purely mechanical description can still capture essential properties of the functional dynamics in proteins. Generally, the current study demonstrates the feasibility of coarse-grained mechanical models in computer experiments of molecular machines. Apparently they cannot replace all-atom simulations. However, given the situation that detailed molecular dynamics methods are not yet capable of describing the entire conformational motions which underlie the transport cycle in ABC transporters, coarser models are potentially useful to explore important aspects of their operation. To study dynamical properties of membrane machines more realistically in computer experiments, it can be attempted in the future to combine elastic-network models with coarse-grained descriptions of the membrane, which would additionally allow to investigate their interplay.

AUTHOR CONTRIBUTIONS
HF designed and performed research, analyzed data and prepared the manuscript.

FUNDING
This work was supported by a Grant-in-Aid for Scientific Research on Innovative Areas 'Spying minority in biological phenomena' (23115007) and Platform for Dynamic Approaches to Living System of the Ministry of Education, Culture, Sports, Science and Technology, Japan. http://www.mext.go.jp/english/. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphy.

2016.00003
Video S1 | Coupled domain dynamics in HmuUV. Conformational motions as predicted by the lowest normal mode (see Materials and Methods). The structure is shown in ribbon representation in front (left) and top view perspective (right).