Multidimensional Umbrella Sampling and Replica-Exchange Molecular Dynamics Simulations for Structure Prediction of Transmembrane Helix Dimers
Structural information of a transmembrane (TM) helix dimer is useful in understanding molecular mechanisms of important biological phenomena such as signal transduction across the cell membrane. Here, we describe an umbrella sampling (US) scheme for predicting the structure of a TM helix dimer in implicit membrane using the interhelical crossing angle and the TM\xe2\x80\x93TM relative rotation angles as the reaction coordinates. This scheme conducts an efficient conformational search on TM\xe2\x80\x93TM contact interfaces, and its robustness is tested by predicting the structures of glycophorin A (GpA) and receptor tyrosine kinase EphA1 (EphA1) TM dimers. The nuclear magnetic resonance (NMR) structures of both proteins correspond to the global free-energy minimum states in their free-energy landscapes. In addition, using the landscape of GpA as a reference, we also examine the protocols of temperature replica-exchange molecular dynamics (REMD) simulations for structure prediction of TM helix dimers in implicit membrane. A wide temperature range in REMD simulations, for example, 250\xe2\x80\x931000 K, is required to efficiently obtain a free-energy landscape consistent with the US simulations. The interhelical crossing angle and the TM\xe2\x80\x93TM relative rotation angles can be used as reaction coordinates in multidimensional US and be good measures for conformational sampling of REMD simulations.
Introduction<!>Definition of the interhelical crossing angle and relative rotation angle<!>Multidimensional US simulations<!>Simulation details<!>Properties of the TM helix dimers of GpA and EphA1<!>US/IMM1 and US/GBSW for structural predictions of GpA and EphA1 TM dimers<!>Comparison between US and REMD simulations for the prediction of GpA structure<!>Conclusions
<p>Membrane proteins, indispensable for many important biological phenomena, function as receptors, transporters, enzymes, and cell recognition and adhesion molecules.[1–5] Although essential for understanding their functions, there are only about 1800 membrane protein structures (400 unique structures) deposited in protein data bank (PDB: http://www.pdb.org) due to experimental difficulties. In particular, the number of atomic structures of transmembrane (TM) helix dimers is less than 20 in the structural databases of membrane proteins like orientations of proteins in membranes (OPM)[6] or protein data bank of transmembrane proteins (PDBTM),[7] although signal transduction across the cell membrane often depends on interactions between TM helix dimers and their conformations.[3,8–10] Therefore, membrane protein structure predictions using computer simulations for membrane protein structure predictions are meaningful if reliable models are constructed within reasonable computational time. In this study, we focus on simulations of TM helix dimers, being the simplest case involving interhelical interactions.</p><p>Glycophorin A (GpA) consists of a single TM helix and forms a tightly packed TM helix dimer in lipid bilayers and micelles.[11,12] GpA has been used as a test system for structure predictions and free-energy calculations of the association/dissociation of a TM helix dimer. As conventional molecular dynamics (MD) simulations are not powerful enough to predict TM helix dimer structures, two enhanced conformational sampling methods have been widely used. They are the replica-exchange method[13–17] (replica-exchange molecular dynamics (REMD) or replica-exchange Monte Carlo (REMC) simulations) and the umbrella sampling (US) method.[18–20] Membrane effects are usually incorporated implicitly using a knowledge-based method[21] or the generalized Born approach[14,22] in these simulations. In addition to these enhanced sampling methods, coarse-grained models have recently been used for simulations of TM helix dimers.[23,24]</p><p>The REMD method has been widely used in the simulations of protein folding, dynamics, and aggregation. In the original method,[13] copies of the target simulation system (replicas) are simulated independently and simultaneously, and every few steps, temperatures and potential energies are exchanged between a pair of replicas. This exchange process induces mixing of temperatures in each replica trajectory and enhances conformational sampling of the system owing to the increased flexibility at high temperatures. But, sometimes the high temperature makes a target molecule too flexible, causing slow convergence of the free-energy calculations. In contrast, the US method[18] or adaptive biasing force simulation[25] requires reaction coordinates for enhanced conformational sampling. For TM helix dimers, interhelical distance, R, is commonly used as a reaction coordinate.[26] This approach successfully provides free-energy changes upon the association/dissociation of a TM helix dimer.[25] However, a single reaction coordinate like R is inadequate for the structural prediction of a TM helix dimer because the contact interfaces between each monomer, which are dependent on the TM-TM crossing angle and TM-TM relative rotation angles (Figure 1), are normally unchanged during simulations.</p><p>Lee and Im, in an effort to solve these problems, developed restrained potentials of the minimum interhelical distance and the interhelical crossing angle for a TM helix dimer.[27] Recently, Park and Im reported that two-dimensional (2D) replica-exchange umbrella sampling (REUS) simulation[28] (this method is also referred to as Hamiltonian REMD[29] or window-exchange umbrella sampling molecular dynamics (USMD)[27,30]) has better performance than temperature REMD or 1D-REUS along R.[27,30] However, the combination of other coordinates like TM–TM relative rotation angles has not been examined so far.</p><p>In this study, we describe a multidimensional US scheme using an interhelical crossing angle (Ω) and TM–TM relative helix rotation angles (ρ) of the TM helix dimer as reaction coordinates. This approach uses the conformation-searching method along the relevant reaction coordinates[31] for TM dimer conformations and has the advantage of MD simulations for free-energy landscapes. To test its applicability and reliability, we carried out US simulations of GpA[11,32] and receptor tyrosine kinase EphA1 (EphA1)[33] TM dimers in implicit membrane. As almost exhaustive sampling on the associated TM helix dimer structure is possible with about 500 windows in US, the free-energy landscapes or the potential of mean forces (PMF) can be used as reference data for designing better protocols for REMD simulations. Therefore, we also performed two sets of REMD simulations in the same implicit membrane model using the GpA TM dimer. A comparison between the PMFs of US and REMD suggests that a wide temperature range in REMD simulations, for example, 250–1000 K, is necessary for efficient structure prediction. Our method is very useful for predicting structures of TM helix dimers without the need of any experimental information other than the amino acid sequences.</p><!><p>Here, we describe an interhelical crossing angle (Ω) and relative rotation angles (ρ) between two helices, both of which are used as the reaction coordinates in the multidimensional US simulations. As described in Figure 1, we define three points in each helix A or B (X1, X2, and X3, where X = A or B) using only the Cα atoms. X1 and X2 are defined as the center of mass (COM) of helix X and the COM of the top half of helix X, respectively. X3 is a selected Cα atom to define the relative rotation angle of helix X with respect to another. Using these six points, we can define the orientation of each helix in a TM helix dimer. The interhelical distance, R, is the distance between A1 and B1. The interhelical crossing angle, Ω, is the dihedral angle of A2–A1–B1–B2: a negative value for the right-handed dimer and vice versa. ρAB is the relative rotation angle for helix A with respect to helix B, and it is defined as the dihedral angle of A3–A2–A1–B1. ρBA is defined in the similar way as ρAB. For GpA and EphA1, these six points are selected as follows: (1)(A1,A2,A3,B1,B2,B3)GPA=(com[X:73−95],com[X:73−84],ca[X:80]),X=A,B)(2)(A1,A2,A3,B1,B2,B3)EphA=(com[X:546−568],com[X:546−557],ca[X:550],X=A,B) where com[X: Y-Z] represents the coordinates of the COM for residues Y-Z in helix X and ca[X: Y] is the Cα coordinate of residue Y in helix X. A3 and B3 are typically chosen within the relatively rigid region or the expected contact motif.</p><!><p>In the multidimensional US method, a two-step procedure is proposed. The first stage is to generate an appropriate initial structure for each US simulation and, in the second stage, the PMF is calculated with the weighted histogram analysis method.[34] First, two ideal helices are placed in parallel to the membrane normal, in the center of the membrane, and separated by 20 Å from each other. Then, the initial orientation of each helix is set to Ω0 and ρ0 at a given window (see below). In the first stage, starting from the separated structure, two TM helices are associated with help of a weak interhelical distance restraint Udis(R): (3)Udis(R)=kR2(R−R0)2 where kR = 2 kcal/mol/Å2 and R0 = 7.5 Å. When the interhelical distance is greater than 11.5 Å2 during the US simulations (in the second stage), a flat bottom harmonic restraint with kR = 0.5 kcal/mol/Å2 is applied to prevent the TM helices from drifting away from each other. The harmonic restraint potential for Ω, Ucross(Ω), and that for ρ, Urot(ρ) are used in both stages: (4)Ucross(Ω)=kΩ2(Ω−Ω0)2,(5)Urot(ρ)=kρ2(ρ−ρ0)2, where kΩ = 500 kcal/mol/radian2 and kρ = 50 kcal/mol/radian2. We used 21-Ω0 values ranging between –50 and 50° in 5° increments, and 24-ρ0 values ranging between 0 and 345° in 15° increments. For a TM heterodimer simulation, when ρAB and ρBA are treated independently, 12,096 (= 21 × 24 × 24) umbrella windows (or independent simulations with different restraint forces) are required. We, therefore, focused only on structure predictions of a homodimer like GpA and EphA1 by assuming the C2 symmetry (i.e., ρ = ρAB = ρBA). In summary, 504 (= 21 × 24) independent short (1 ns) simulations with different restraints for Ω and ρ in the first stage and long (40 ns) simulations in the second stage were carried out to obtain the 2D PMF map along Ω and ρ.</p><!><p>All the simulations were performed with the CHARMM package[35] and the multiscale modeling tools for structural biology (MMTSB) toolset.[36] We used two implicit solvent/membrane models, IMM1[37,38] and generalized-born with a simple switching function (GBSW),[14] and utilized the united-atom CHARMM C19 EEF1.1 force field[39,40] for IMM1 and all-atom CHARMM C22 force field[41] with dihedral cross-term corrections (CMAP)[42] for GBSW. The IMM1 model is based on the Gaussian solvent-exclusion approximation and takes 10 Å of solvent/membrane interface between the interior cyclohexane region and the exterior water region and a 27 Å width of interior region to mimic a POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) membrane bilayer. In the GBSW model, we used 0.03 kcal/mol/Å2 for surface tension, 5 Å for each membrane/solvent interface, and 27.5 Å for the thickness of the membrane hydrophobic region to mimic the same membrane The nonbonded interactions were switched off over 7–9 Å for IMM1, whereas a 16 Å cutoff distance for nonbonding interactions was applied for GBSW. In both models, the planar membrane is perpendicular to the Z-axis and centered at Z = 0. Langevin thermostat[43,44] controlled the temperature at 300 K. SHAKE[45] was applied to the covalent bonds involving hydrogen atoms with a time step of 2 fs. Hereafter, the multidimensional US simulations with IMM1 and GBSW membrane models are denoted as US/IMM1 and US/GBSW, respectively.</p><p>We also performed REMD simulations of the GpA dimer with IMM1 implicit membrane. In one simulation, 32 replicas were distributed exponentially over a temperature range of 250–1000 K (REMD250–1000K), whereas another consisted of 16 replicas distributed exponentially over a temperature range of 300–550 K (REMD300–550K). Replica-exchange between a pair of replicas was carried out every 2 ps. We simulated for 100 ns for each replica in REMD250–100 K and REMD300–550 K, respectively. To avoid artificial drifts of the GpA helix dimer toward solvent regions and its unfolding at high temperatures, we added a flat bottom harmonic restraint potential (k = 10 kcal/mol/Å2) to the COM of each helix when it is 5 Å from Z = 0 and Cα-RMSD harmonic potentials (10 kcal/mol/Å2) for each helix (as observed in the PDB structures). In addition, we added a flat bottom harmonic restraint potential (kR = 0.5 kcal/mol/Å2) to the GpA dimer when R > 11.5 Å. Other simulation parameters were the same as US/IMM1.</p><!><p>We first summarize the structural properties of GpA[11,32] and EphA1[33] TM helix dimer structures in PDB. In Figure 2A, the amino sequences of GpA and EphA1 used in this study are shown. There are two "G(A)xxxG(A)" sequence motifs in GpA (G79xxxG83, A82xxxG86) and four in EphA1 (G546xxxA550, A550xxxG554, G554xxxG558, A560xxxG564). The "G(A)xxxG(A)" motif is often found at TM–TM contact interfaces.[10,46–48] In addition, β-branched residues (Val, Leu, Ile) adjacent to the GxxxG-like motif (called GVxxGV-like motif) enhance the stability of the GpA helix dimer.[47] Averaged structures derived from 60 NMR structures for GpA (1AFO,[11] 2KPE,[32] and 2KPF[32]) and 12 NMR structures for EphA1 (2K1K[33]) are shown in Figure 2B. The contact residues are defined as those within 6 Å between Cα–Cα atoms. In the average structure of GpA, the contact residues are L75, I76, G79, V80, G83, V84, and T87, and in that of EphA1, they are G546, E547, A550, V551, G554, L555, G558, and A559. The contact motifs "LIxxGVxxGVxxT" for GpA and "AxxxGxxxG" for EphA1 fall within the defined distance limits. Gly and Ala residues located at the contact interfaces are shown in red in the sequence in Figure 2A. Note that "A82xxxG86" in GpA and "A560xxxG564" in EphA1 are not at the contact interfaces. The interhelical coordinates (R, Ω, ρAB, ρBA) calculated from the GpA and EphA1 NMR structures are listed in Table 1. There are large apparent deviations for three PDB structures of GpA presumably due to the different experimental environments of the TM helix dimers.[11,32,33]</p><!><p>Figure 3 shows the 2D-PMFs for the GpA TM dimer from US/IMM1 and US/GBSW as a function of interhelical crossing angle (Ω) and relative helix rotation angle with respect to one another (ρ). The length of the trajectory used in the PMF calculation in the two membrane models is estimated from the convergence of the free energy difference between the left-handed and right-handed conformations shown in Supporting Information Figure S1. 1–40-ns simulation trajectories were used to calculate the free energy surface for IMM1, whereas 14–40-ns simulation trajectories served in the PMF calculations for GBSW. The lowest free energy is located at (Ω, ρ) = (–45,–16) and (–40,–18) for IMM1 and GBSW, respectively. Both structures adopt right-handed conformations. We estimated block errors (standard deviation) for the left-handed dimer for the block length of 1 ns, and the error was 0.3 kcal/mol in both models. The average structure of the lowest free energy for the US window for the two implicit membrane models is shown on the left in Figure 3; the structures have Cα RMSD of 1.8 Å (US/IMM1) and 1.5 Å (US/GBSW) from the average NMR structure, respectively. The contact residues are indicated in Figure 3. The right-handed dimer is the most stable conformation, and L75, I76, G79, V80, G83, V84, and T87 are located at the contact interface and form a "LIxxGVxxGVxxT" contact motif as in the NMR structures.[11,32]</p><p>The 2D-PMF maps for EphA1 from US/IMM1 and US/GBSW are shown in Figure 4. In US/IMM1, the helix tended to partially unfold in some windows, so an RMSD restraint was applied to the Cα atoms in each helix. The convergence of free energy differences between the left-handed and right-handed conformations is shown in Supporting Information Figure S2. Accordingly, 17–40-ns simulation trajectories were used to calculate the free energy surface for US/IMM1, whereas 27–40-ns simulation trajectories were used to calculate that for US/GBSW. The average structure with the lowest free energy adopts a right-handed conformation (shown on the left-handed side in Fig. 4). The Cα RMSDs between the simulated structures and the NMR structure are 0.8 Å (US/IMM1) and 1.1 Å (US/GBSW), respectively. Right- and left-handed conformations have similar free energies (0.59 and 0.37 kcal/mol differences for US/IMM1 and US/GBSW, respectively). The contact residues are A550, V551, G554, L555, and G558, including the glycine zipper motif AxxxGxxxG, as also seen experimentally by NMR.[33] Left-handed dimer structures are shown on the right in Figure 4, and their interhelical crossing angles are small. The contact interface of left-handed dimers contains either a "G546xxxA550" (a GxxxG-like) motif or a "L556xxxxxxL563" leucine zipper motif. We suspect that there is a second possible left-handed dimer conformation with the leucine zipper contact motif. EphA2, another EphA TM dimer, adopts a left-handed conformation and uses the leucine zipper, instead of the glycine zipper motif, as a contact motif.[49] Note that if a dimer sequence contains many wellknown contact motifs such as in EphA1, the convergence of PMF is much slower than for the simpler case of GpA.</p><p>For both TM dimers, the NMR structures occur at the lowest free-energy states of the PMF landscapes. The differences in barrier height and location are related to the balance between electrostatic interactions and solvation energy,[50] in which case they would originate from the plasticity of the system due to variation in the force field parameters of proteins.[39–42]</p><!><p>In the previous simulations of GpA using the temperature REMD method in implicit membranes, it was difficult to obtain an NMR structure with the right-handed conformation as the global minimum.[14] Im et al. have successfully simulated the NMR structure of GpA using the two-fold symmetry.[14] The free-energy landscapes of GpA obtained with our simulation schemes (2D-US and REMD) allow us to verify our protocols for REMD simulations with IMM1. We chose the temperature ranges 250–1000 K (REMD250–1000K) and 300–550 K (REMD300–550K) for the REMD simulations. A harmonic RMSD restraint potentials to Cα atoms in each helix was used to prevent unfolding of each TM helix at the highest temperatures. This approach is a similar to that used in previous REMC simulations of rigid-body translational movements of TM helices.[15]</p><p>First, random walks in replica, temperature, and potential energy spaces for REMD250–1000K (Fig. 5) and REMD300–550K (Supporting Information Fig. S3) show that the REMD algorithm works well in both simulations. Figure 5 and Supporting Information Figure S3 also show random walks in conformational spaces, like RMSD from the averaged NMR structure, the interhelical distance (R), interhelical crossing angle (Ω), and relative rotation angles (ρAB and ρBA). In the REMD250–1000 K trajectory, random walks are observed in all structural properties, and at 40 and 90 ns, they correspond to those of the NMR structure. In other time frames, they take different values, suggesting that the free-energy landscape has a funnel feature similar to soluble globular proteins.[51,52] In contrast, in REMD300–550K random walks in relative rotational angles are slower than in REMD250–1000K. Although ρAB and ρBA in REMD300–550K show values ranging widely, higher temperatures seem to be required for better conformational sampling.</p><p>We also calculated the 2D-PMF along with Ω and ρ from the trajectories of REMD simulations. Note that we did not use the symmetry restraint for each TM helix in both REMD simulations and we only plot ρAB of GpA dimer in Figure 6. In Figure 6A, 2D-PMF at 300 K in REMD250–1000 K has several common features with that in US/IMM1 or US/GBSW (see Fig. 3). First, we observe the global free-energy minimum around (Ω, ρ) = (–38, –18) in Figure 6A, where both ρAB and ρBA of each TM helix take almost the same values. Thus, the right-handed NMR structure of GpA TM helix appears as the most stable conformation in REMD250–1000K. In Figure 6E, we observed four representative structural clusters calculated from k-means clustering implemented in the MMTSB toolset.[36] The representative conformation is mapped to the 2D-PMF at 300 K in REMD250–1000K. The second stable conformation in US, which is the left-handed conformation (LH1, in Fig. 3A), corresponds to the cluster 2 and 3. We also compare the 2D-PMF at the high temperatures in REMD250–1000K (550, 750, and 1000 K) in Figures 6B–6D. The large interhelical crossing angle for GpA is available only at lower temperature due to molecular interactions at the contant interface. We also show the 2D-PMF at 300 and 550 K in REMD300–550K in Supporting Information Figure S4. Although a wider conformational sampling was achieved at 550 K, the contact interfaces, which are determined with ρ, are correlated with the crossing angle, suggesting that 550 K may not be sufficiently high for taking all the possible contact interfaces. The relative free-energy values and clustering results in REMD300–550K are slightly different from those in REMD250–1000K.</p><p>Finally, we discuss appropriate protocols of multidimensional US and temperature REMD simulations for TM helix dimers in implicit membranes. Multidimensional US used with the interhelical crossing angle and TM–TM relative rotation angles as reaction coordinates is a promising tool for predicting "homo"-dimer structures. All the contact interfaces for a TM helix dimer can be analyzed. As indicated in the Method section, about 500 umbrella windows need to be prepared for the exhaustive searches of associated structures. Such simulations can be done with a massively parallel supercomputer. Application to a TM "hetero"-dimer requires a much greater number of windows (12,096). Computing such numbers is usually out of reach using current computational resources so that additional methodological developments to reduce the number of windows are required, such as replica-exchange metadynamics.[53] It is noted that the current method is not developed to measure the TM–TM interaction energy, and such calculations require the PMF calculations along the TM–TM distance. However, one may expect a flat landscape instead of a well-defined funnel-like landscape on the 2D-PMF map using the current method for a noninteracting dimer. Our results emphasize the importance of a wider range of temperatures for successful REMD simulations. The PMF landscapes along Ω and ρ are useful for judging whether there is sufficient conformational sampling with REMD. This landscape is also beneficial for judging long MD simulations of a TM helix dimer using explicit solvent and membrane models[54,55] as well as coarse-grained models.[24,56] The computational cost of REMD is much less than that required in the multidimensional US. The total simulation time in the multidimensional US for GpA or EphA1 helix dimers is 20.664 μsec = (1 + 40) ns × 504 windows, whereas in REMD, it is just 3.2 μsec 5 100 ns × 32 replicas for GpA. This suggests that REMD may be still useful for structure prediction of a TM helix dimer, provided that appropriate simulation protocols are selected. The Cα-RMSD restraints in this study were necessary to keep the helical structures (as observed in the PDB structures) during the simulations. This issue is more related to the development of the proper force field to capture such conformations. Having a good sampling method such as the one we developed here could in principle help such developments. However, a weaker Cα-RMSD restraint can be applied in the current method in the consideration of helix flexibility.</p><!><p>We have described a multidimensional US scheme for the prediction of TM helix dimer structure in biological membranes. It uses the interhelical crossing angle and TM–TM relative interhelical rotation angles as reaction coordinates to search different contact interfaces between TM helix dimers. By assuming that TM-homo dimers have symmetrical structures, we can reduce the number of US windows into about 500. We tested this approach to predict the structures of GpA and EphA1 TM dimers in two different implicit membrane models (IMM1 and GBSW). The simulated structures at the lowest free-energy minima correspond to the averaged NMR structures and agree with the results of REMD simulations at 250–1000 K. Conformational sampling along the interhelical crossing angle and TM–TM relative rotational angles is found to be essential for successful REMD simulations of TM helix dimers. These angles are also useful to examine the convergence of long MD simulations of a TM helix dimer in explicit solvent and lipid-bilayer or in coarse-grained models.</p>