Behavior of Solvent-Exposed Hydrophobic Groove in the Anti-Apoptotic Bcl-XL Protein: Clues for Its Ability to Bind Diverse BH3 Ligands from MD Simulations

Bcl-XL is a member of Bcl-2 family of proteins involved in the regulation of intrinsic pathway of apoptosis. Its overexpression in many human cancers makes it an important target for anti-cancer drugs. Bcl-XL interacts with the BH3 domain of several pro-apoptotic Bcl-2 partners. This helical bundle protein has a pronounced hydrophobic groove which acts as a binding region for the BH3 domains. Eight independent molecular dynamics simulations of the apo/holo forms of Bcl-XL were carried out to investigate the behavior of solvent-exposed hydrophobic groove. The simulations used either a twin-range cut-off or particle mesh Ewald (PME) scheme to treat long-range interactions. Destabilization of the BH3 domain-containing helix H2 was observed in all four twin-range cut-off simulations. Most of the other major helices remained stable. The unwinding of H2 can be related to the ability of Bcl-XL to bind diverse BH3 ligands. The loss of helical character can also be linked to the formation of homo- or hetero-dimers in Bcl-2 proteins. Several experimental studies have suggested that exposure of BH3 domain is a crucial event before they form dimers. Thus unwinding of H2 seems to be functionally very important. The four PME simulations, however, revealed a stable helix H2. It is possible that the H2 unfolding might occur in PME simulations at longer time scales. Hydrophobic residues in the hydrophobic groove are involved in stable interactions among themselves. The solvent accessible surface areas of bulky hydrophobic residues in the groove are significantly buried by the loop LB connecting the helix H2 and subsequent helix. These observations help to understand how the hydrophobic patch in Bcl-XL remains stable in the solvent-exposed state. We suggest that both the destabilization of helix H2 and the conformational heterogeneity of loop LB are important factors for binding of diverse ligands in the hydrophobic groove of Bcl-XL.


Introduction
The role of Bcl-2 family in regulating the mitochondrial outer membrane permeabilization, an important step in the cell death process, is well established [1,2,3,4,5]. Among the functionally apposing Bcl-2 members, the anti-apoptotic Bcl-X L protein is one of the first proteins from the Bcl-2 family that has been investigated by several researchers [6,7,8,9,10]. This protein has been specifically shown to be over-expressed in several human cancers including lung [11], colon [12], pancreatic [13], breast [14] and prostate [15] cancers. Hence, Bcl-X L like other antiapoptotic Bcl-2 members is an attractive target for anti-cancer drugs [16,17]. Structural knowledge of this protein has helped to understand how Bcl-X L is recognized by its pro-apoptotic partners which is a crucial step for the cell survival [18]. The pronounced hydrophobic groove present in this helical bundle protein acts as a binding region for the BH3 domain of pro-apoptotic Bcl-2 members. Experimental studies have shown that the BH3 peptides derived from the pro-apoptotic proteins have the ability to elicit the same biological response as that of the parent proteins [19]. As of December 2012, 50 Bcl-X L structures have been deposited in the Protein Data Bank [20] including 30 structures in which Bcl-X L is in complex with a BH3 peptide or an inhibitor molecule. However, biochemical and pharmacological studies have demonstrated that different pro-apoptotic proteins exhibit distinct binding affinities for Bcl-X L although experimentally determined structures show that they all bind to the same hydrophobic groove [19,21,22]. Structure-based drug design approach combined with docking simulations has been used to identify and/or design small molecule inhibitors that can target Bcl-X L [23]. Design of molecules mimicking the BH3 domains of pro-apoptotic Bcl-2 proteins [24,25,26], pharmacophore-based database searching [27,28] and fragment-based building of bioactive inhibitors [29] are some of the approaches used in designing and synthesizing Bcl-X L -specific inhibitors.
Three factors could be important for the pro-apoptotic proteins to recognize and bind particular anti-apoptotic Bcl-2 proteins with specific binding affinities.
(i) Specific interactions could help the ligands to bind tightly to the target proteins. Among the different computational approaches that investigated the differential binding affinities, molecular dynamics (MD) simulations have identified certain hydrophilic interactions between residues from Bcl-X L protein and BH3 peptides [30]. These interactions involve the Arg residue next to the highly conserved Leu residue in the BH3 peptide and the Glu residue present in the helix preceding the central hydrophobic helix of Bcl-X L protein. Such interactions are suggested to be responsible for the high binding affinity of BH3 peptides of certain proapoptotic proteins and the results are supported by affinity studies.
(ii) Spectroscopic and mutational studies have demonstrated a possible correlation between the BH3 peptides and their helical stability [31]. MD simulations of several Bcl-X L -and A1-binding BH3 peptides showed that the BH3 peptides that have high-binding affinities for either of the protein exhibited stable helical segments in aqueous medium [32,33]. Analysis of the simulated structures revealed the reasons for the stable nature of the amphipathic BH3 peptide helices. Although clustering of hydrophobic residues destabilized the helices to some extent, capping interactions and stable hydrophilic interactions from the solvent-exposed hydrophilic side of the amphipathic BH3 helices helped to retain the helical nature of high-affinity BH3 peptides in the solvent. Such studies can aid in designing BH3 mutants and BH3 mimetics that can bind to the anti-apoptotic proteins with very high binding affinity. (iii) The last and the most important factor is the flexibility in the hydrophobic-binding region. Since Bcl-X L is recognized by the BH3 region of several pro-apoptotic proteins [5,19,22,31,34], it is important that the hydrophobic groove must accommodate diverse BH3 sequences with some common features like conserved hydrophobic residues. Very few studies have looked into this aspect of the Bcl-2 proteins. Moreover, it is intriguing to note that such a large hydrophobic patch is exposed when the BH3 peptide is not bound to the protein. To our knowledge, the behavior of solvent exposed hydrophobic groove is not thoroughly explored. In this study, we have carried out molecular dynamics simulations of Bcl-X L protein in the apo-form. We have also considered structures of Bcl-X L complexes in which the bound BH3 peptides were removed (holo Bcl-X L ). Our results from multiple MD simulations of Bcl-X L in apoand holo-forms show that the hydrophobic residues from Bcl-X L which were interacting with the bound BH3 peptides in the structures of complexes interact among themselves in apo and holo structures. Accessible surface area calculations indicate that the exposed bulky hydrophobic residues in the hydrophobic groove are partially shielded by the loop connecting the BH3-containing helix and the next helix. Our simulations using twin-range cut-off exhibit destabilization of BH3-containing helix whereas simulations which used Particle Mesh Ewald (PME) have shown the same helix as stable. Experiments have indicated that the BH3containing helix has to undergo functionally important conformational changes for homo-or hetero-dimerization [35,36]. Unwinding of BH3-containing helix can also help in binding diverse BH3 ligands derived from pro-apoptotic Bcl-2 proteins. It is possible that the event of BH3containing helix destabilization in PME simulations requires a longer time-scale (from several hundred nanoseconds to microseconds time scale).

Initial structures
The starting structure for apo-Bcl-X L corresponds to the PDB ID 1PQ0 (Resolution: 2.2 Å´) [37]. The initial structure of holo-Bcl-X L was obtained from the structure of Bcl-X L /Bim complex (PDB ID: 1PQ1; resolution: 1.65 Å´) by removing the bound Bim peptide. In holo structure, by removing the bound BH3 peptide, the previously buried hydrophobic groove is now exposed. MD simulations of holo structure will give an idea whether the behavior of exposed hydrophobic patch is similar to that observed in apo structure simulations. Moreover, one will be curious to know what happens to those hydrophobic residues which were originally interacting with the BH3 peptide ligand in the structure of the complex. Hence simulations were carried out on both apoand holo-Bcl-X L structures and the results are compared.
It must be noted that in apo and holo structures, the disordered long loop linking the first two helices is not fully resolved and the structures were determined by deleting the C-terminal transmembrane domain [37]. The structure of the missing loop was built using the ''Homology'' module of InsightII suite of software (Accelrys Inc., San Diego, CA) by considering another structure of Bcl-X L complex (PDB ID: 1G5J [31]) as template. Wherever sidechain atoms are not resolved in the experimental structures, they were constructed using the ''Biopolymer'' module of InsightII. Residue numbers of Bcl-X L correspond to the sequence with UniProt [38] accession code Q64373. The ribbon diagram of a Bcl-X L structure with its different helices is shown in Figure 1.
In the PDB structures, hydrogen atoms were built using the utility tool PDB2GMX available in the GROMACS modeling software version 3.2.1 [39,40]. Only those hydrogens bonded to the polar atoms or aromatic carbon atoms were added. The Nand C-terminal residues were capped using the -NH 2 and -COOH groups respectively. All titratable residues were assumed to have default protonation states. The starting structures were placed at the center of a cubic box whose size was determined such that the distance between the protein and the edge of the box was Figure 1. Ribbon representation of Bcl-X L helical bundle structure. The major helices are labeled as H1 to H6 in this apostructure (PDB ID: 1PQ0) and are shown in different colors. The loops (LB, LC and LD) connecting different helical segments are shown in purple color. The loop LA connecting helices H1 and H2 is not shown for clarity. Molecular plots in this figure and subsequent figures were generated using the VMD software [84]. doi:10.1371/journal.pone.0054397.g001 Dynamics of Hydrophobic Groove in Bcl-X L PLOS ONE | www.plosone.org at least 10 Å´. SPC [41] water molecules were used to solvate the Bcl-X L protein using GENBOX tool in GROMACS. Energy minimization of solvated protein was carried out using steepest descent and conjugate gradient methods. Harmonic positional restraints were applied on protein heavy atoms during minimization.

Simulation protocol
Eight independent simulations of the solvated and energyminimized Bcl-X L structures were carried out for a period of 25 to 55 ns using GROMOS96 43 a1 force field [42,43]. The simulations differed in the starting structures (apo-or holo-Bcl-X L ; see above), equilibration protocol (two different schemes) and in the treatment of long-range interactions (twin-range cut-off or PME methods). In the simulations, the bonds and angles of solvent molecules were constrained using SETTLE algorithm [44]. All the bonds involving the hydrogen atoms of the solute were constrained using LINCS algorithm [45]. Periodic boundary conditions were applied in all three directions. For all the simulations, the reference temperature and pressure were considered to be 300 K and 1 bar respectively. Simulation temperature and pressure were maintained using Berendsen's algorithm [46]. For temperature, a coupling constant of 0.01 ps during equilibration and 0.05 ps for the production run was used. The boundary pressure was maintained with a coupling constant of 1.0 ps. The solute and solvent temperatures were coupled separately. A time-step of 0.002 ps was used and the non-bonded list was updated every 10 steps.

Two different treatments of long-range interactions
In simulations that used twin-range spherical cut-off, nonbonded interactions were evaluated using a cut-off of 10 and 18 Å´. In PME simulations, the real space electrostatic and van der Waals interactions were evaluated with a spherical cut-off of 10 Å´. A grid space of 1.5 Å´was used with sixth order interpolation. The system was neutralized by adding 13 Na+ ions at random positions by replacing water molecules. The GENION tool in GROMACS was used for this purpose.

Equilibration Scheme I
The system was first equilibrated for 1 ns with harmonic positional restraints applied on the heavy atoms of Bcl-X L protein using NVT ensemble (constant number of atoms, volume and temperature). In the next stage, all restraints were removed and the system was equilibrated for another 1 ns using NPT (constant number of atoms, pressure and temperature) ensemble. After 2 ns of equilibration, the apo-and holo-Bcl-X L structures were simulated for a period of 25 to 55 ns using the same NPT ensemble.

Equilibration Scheme II
In the simulation of Bcl-X L complexes, we have previously shown unwinding of helix H2 [30]. A similar observation was made in simulations in which the systems were equilibrated using Scheme I (see Results section). One can argue that such unwinding could be due to artifacts resulting from inadequate equilibration of the systems. Hence in this scheme, the equilibration protocol was modified. After the first 1 ns of equilibration in which harmonic positional restraints were applied on all heavy atoms of the protein, the next 1 ns equilibration was carried out in four stages with restraints imposed on helix H2 and its environment. In the first stage, 250 ps of equilibration run was carried out by imposing harmonic positional restraints on residue pairs from helix H2 and other helices (H1, H4 and H6) if these residues have at least one pair of heavy atoms with distance within 4 Å´. In the next three stages each consisting of 250 ps equilibration run, positional restraints were imposed first on all heavy atoms of H2, then only on the backbone atoms of H2 and finally only on the Ca atoms of H2. The last 1 ns of equilibration was carried out without any restraints. Out of a total of 3 ns equilibration, the first 2 ns was performed using NVT ensemble. The last 1 ns of equilibration and the production run used NPT ensemble. The apo-and holostructures were simulated for a period of 55 ns production run after the 3 ns equilibration. All the apo-and holo-Bcl-X L simulations using two different equilibration schemes and two different treatments of long-range interactions are summarized in Table 1.

Results
We have presented and compared the results of all eight simulations for both apo-and holo-Bcl-X L systems. Average helical content of Bcl-X L , stability of individual helical segments, root mean square deviation (RMSD) with respect to the starting structure, solvent accessibility and interactions of hydrophobic residues in the hydrophobic groove are some of the analyses that have been carried out for different simulations.

Stability of Bcl-X L protein
Helical content of Bcl-X L . The average helical content of Bcl-X L was calculated for each MD simulated structure as described in our earlier studies [30]. First, Q, y values were calculated for each residue and then a residue is considered to adopt right-handed a-helical conformation if the backbone dihedral angles satisfy the following criterion [47] {140 0 ƒqƒ{30 0 and {90 0 ƒyƒz45 0 ð1Þ Average helical content was calculated as the fraction of residues in Bcl-X L having helical conformation. MD trajectories of average helical content for the four twin-range cut-off simulations are plotted in Figure 2A as the function of time. It is clear that the helical content of Bcl-X L decreases to about 65 to 70% compared to almost 80% observed in the experimentally determined structures. The helical content is stabilized during the last 20 ns of production run. Compared to apo and holo simulations, Bcl-X L retained a larger fraction of helical structure (75%) when it is bound to BH3 peptide of pro-apoptotic proteins as evident from the simulations of Bcl-X L complexes [30]. Helix content of each MD simulated structure was calculated and the MD trajectories of all four PME simulations are displayed in Figure 3A. It is clear that the helix content is maintained close to the experimentally determined structure (,80%) in all the PME apo-and holo-Bcl-X L simulations. This is in contrast to what was observed in twin-range cut-off simulations in which the protein lost about 10 to 15% of the total helical content in both apo and holo simulations. This demonstrates that the Bcl-X L helix bundle structure is more stable in PME simulations compared to the twinrange cut-off simulations for the time period of simulation.

Stability of individual helical segments
The decrease in the helical content of Bcl-X L in twin-range cutoff simulations indicates that some of the helices might have lost the helical conformation either partially or fully. To identify and analyze the individual helical segments in the protein, we used the following definition. A segment in the polypeptide chain Bcl-X L was considered helical if Q and y values of at least four consecutive residues must be in helical conformation. In other words, the Q and y values of at least four consecutive residues must satisfy equation (1). This definition was applied to identify the helical segments in the experimentally determined and MD simulated structures. We also defined a helical segment as stable if at least four consecutive residues maintain helical conformation for at least 80% of the simulation time.
Our criteria identified six helical segments designated as H1 to H6 in the experimentally determined structures ( Table 2). The loop LB connecting H2 and H3 also contains a small helical segment. We also identified stable helical regions in the  Figure 2. Stability of Bcl-X L structure in twin-range cut-off simulations. MD trajectories of (A) percentage helical content and (B) RMSD profiles shown for the four simulations that used twin-range cutoff to evaluate the long-range interactions. Percentage helical content of each MD simulated structure was calculated as described in the Methods section. RMSD was calculated by considering the Ca atoms of all stable helical segments H1 to H6 (Table 2). doi:10.1371/journal.pone.0054397.g002 Figure 3. Stability of Bcl-X L structure in PME simulations. MD trajectories of (A) percentage helical content and (B) RMSD profiles shown for the four simulations that used PME to treat the long-range interactions. Percentage helical content of each MD simulated structure was calculated as described in the Methods section. RMSD was calculated by considering the Ca atoms of all stable helical segments H1 to H6 ( simulations Apo-I, Apo-II, Holo-I and Holo-II. We identified helical segments which are stable in at least two out of four twinrange cut-off simulations ( Table 2). We noticed the formation of a small stable helical segment of just one turn in the loop LA connecting the helices H1 and H2. It is evident that there is a loss in helicity in helix H2 and this helix is broken in the middle ( Figure  S1). In fact in two simulations (Apo-I and Holo-II), this helix is almost completely destabilized. A similar observation was made in the simulations of three Bcl-X L complexes although the unwinding was less dramatic in those cases [30]. This unwinding was also shown in a simulation using a different force field [48] with spherical cut-off [30]. We then examined the stability of helix H2 in all PME simulations and not surprisingly, this helix was found to be very stable (See Figure S1). This observation is unlike what was observed in the twin-range simulations in which significant portion of helix H2 lost the helical character.
As far as the other helices are concerned, with the exception of the C-terminal helix H6, other major helices are more or less intact throughout all the four twin-range cut-off simulations. The terminal regions of a polypeptide chain are usually flexible and in this case, the region after helix H6 is truncated. Hence, the instability in helix H6 was not surprising.
There appears to be several reasons for the destabilization of helix H2 in twin-range cut-off simulations. In the complex simulations, when a similar observation was made, it was suggested that the destabilization of helix H2 could be due to the presence of a glycine residue (Gly 94), or presence of residue pairs which are either acidic [(E92, D95), (E92, E96), (D95, E98)] or basic (K87, R91). These residue pairs are separated by three to four residues [30]. In such an arrangement, the charged residues will be one above the other if they are present in an a-helix and as a result, these like-charged residue pairs will repel each other. A third factor was suggested to be the absence of stable inter-helix interactions involving helix H2 [30]. In Bcl-X L apo/holo simulations as in complex simulations, we have analyzed all three factors (data not shown). Based on this analysis, we reached a similar conclusion in apo/holo Bcl-X L simulations that all three factors individually or together could contribute to the destabilization of helix H2. MD simulated structures saved at the end of the production runs superimposed on the starting structure for the apo (Apo-I and Apo-II) and holo (Holo-I and Holo-II) simulations are shown in Figure 4A and 4B respectively. In PME simulations, the length of the simulations perhaps could be the reason that the unwinding has not been observed in helix H2. Longer simulations could have resulted in destabilization of the same helix. MD simulated structures superposed with the starting structure are shown in Figure 5 for the PME simulations.

RMSD analysis
MD trajectories of RMSD (Root Mean Squared Deviation) calculated between the simulated structures and the starting structure are plotted in Figure 2B for the four twin-range cut-off simulations. Only C a atoms of the residues from the six stable helical segments, H1 to H6 (Table 2), were considered for this purpose. At the end of 55 ns simulations, the structures deviated from the starting structures by close to 5 Å´. This deviation is much larger, at least by 1.5 to 2 Å´, compared to the Bcl-X L complex simulations [30]. This indicates that the interacting BH3 peptides restrict the overall dynamics of the Bcl-X L protein. In both apo and holo forms, Bcl-X L does not have any constraints in the form of intermolecular interactions and the protein is relatively more flexible as evident from the RMSD analysis.
We next calculated RMSD of each MD simulated structure with respect to the starting structure for the PME simulations. As in the twin-range simulations, only the Ca atoms of the six major helical segments (H1 to H6; see Table 2) were considered for this purpose. Analysis of RMSD trajectories shows that all PME simulations exhibit a very low RMSD of 1.5 Å´( Figure 3B). This again demonstrates that the protein has deviated less from the starting structures in all PME simulations.
Interactions of solvent-exposed hydrophobic residues in the hydrophobic groove While the hydrophobic groove is shielded from the solvent by the bound BH3 peptide in structures of complexes, it is exposed in apo-Bcl-X L . In general, exposure of a large hydrophobic patch in aqueous medium is not energetically favorable. Hence, it is   intriguing to find out the behavior of this hydrophobic patch in the solvent-exposed state in both apo-and holo-Bcl-X L structures. First, we specifically focused our attention on those hydrophobic residues which participated in interactions with the bound BH3 peptide in structures of complexes. We analyzed whether these residues now participate in interactions with other residues. For this purpose, we identified 174 interacting residue pairs from four experimentally determined structures; one apo-Bcl-X L (PDB ID: 1PQ0) and three holo-Bcl-X L structures (PDB IDs: 1BXL, 1G5J and 1PQ1; the bound BH3 peptides were removed from the structures of complexes). We define two residues as interacting residues if at least one pair of heavy atoms between the residues is within 4 Å´. The interacting residue pairs occur between the six helices H1 to H6 as inter-helical interactions or between one of these helices and loops LB, LC or LD (Figure 1). We included these loops since they also contribute to the hydrophobic groove. The minimum distances between these residue pairs were followed throughout the simulations and analyzed for the last 20 ns of the production runs. An interaction is defined as stable if the minimum distance between at least one pair of heavy atoms is less than 4 Å´for more than 50% of the analysis period in at least 2 out of 4 twin-range cut-off simulations. We have identified a total of 38 residue pairs which can be considered as participating in stable interactions. Majority of 38 interactions are observed in two out of three Bcl-X L complex simulations [30] and are common with apo/holo simulations (data not shown). These interactions are mostly hydrophobic and occur in the interior of the helix bundle providing overall structural stability for the protein.
Eight stable interactions were identified only in the apo-and holo-Bcl-X L simulations ( Table 3). Five of them involve residue pairs in which at least one residue participated in stable interactions with a BH3 peptide in structures of complexes [30]. The residues from these residue pairs are from helices and loops that form the hydrophobic groove. This shows that in the absence of a bound BH3 peptide ligand, the exposed hydrophobic residue from the hydrophobic groove can compensate the energy penalty to some extent by interacting with another hydrophobic residue in the same hydrophobic groove. The remaining three stable interactions involve a residue in loop LB which is known to interact with the BH3 peptide [30]. All the eight stable interactions are shown in Figure 6A and 6B from Apo-I and Holo-I simulations. The same data for the other six simulations are provided in Figures S2, S3, S4. It is clear that these stable hydrophobic interactions which are absent in the Bcl-X L complex simulations are spread throughout the hydrophobic groove from one end to the other end.

Solvent accessibility of hydrophobic residues in the hydrophobic groove
We calculated the solvent accessible surface area (SASA) of sixteen hydrophobic residues (A89, A93, F97, V126, L130, V135, V141, A142, F146, L150, M159, L162, W181, F191, L194 and Y195) that participate in the formation of hydrophobic groove and majority of these residues were in stable contact with at least one BH3 peptides (Bak, Bad and Bim) identified in our earlier simulations studies [30] of Bcl-X L complexes. In the apo-and holo-Bcl-X L structures, these residues are free, exposed and there is no bound ligand to interact. Since loop LB wraps around the binding BH3 peptides or ligands in the structures of complexes, SASA values for these hydrophobic residues were calculated with and without the loop LB (residues 101 to 119). This will give an idea of how much the loop LB contributes in burying the exposed hydrophobic residues. We also performed the same analysis for all the six hydrophobic residues (Y101, A104, F105, L108, L112 and I114) from the loop LB. The GROMACS (version 3.3) utility tool g_sas was used for this purpose. Average SASA values were calculated for the last 5 ns of the production run (Table 4) for the residues forming the hydrophobic groove. For loop LB hydrophobic residues, average SASA values have been reported for the first and the last 5 ns of the production run (Table 5).
We found that average SASA values of 8 out of 16 residues (A89, A93, V135, V141, W181, F191, L194 and Y195) are not affected by loop LB (data not shown). These residues are present in the unwound region of helix H2 or close to other loop regions. Accessibility of the seven residues (F97, V126, L130, F146, L150, M159 and L162) is significantly affected by the presence of loop LB. These are bulky hydrophobic residues and the loop LB buries up to 75 Å´2 surface area for each of these seven residues. They are present in the regions (helices H2, H3 and H4 and loop LD connecting helices H4 and H5) that significantly contribute to the formation of hydrophobic groove ( Figure 6C and 6D). It is also interesting to note that six out of seven of these residues also participate in stable hydrophobic interactions only in apo-and holo-Bcl-X L simulations (Table 3). In all four simulations, presence of loop LB helps to bury a total of 221 to 380 Å´2 in these seven hydrophobic residues ( Figure 6E and 6F; also see Figures S2, S3, S4).
We next analyzed the SASA values of all the hydrophobic residues from loop LB to see whether a specific pattern is observed in burying or exposing the hydrophobic residues across the eight simulations. In all the twin-range cut-off simulations, the Y101 residue which was initially exposed to the solvent was buried at the end of 55 ns production run (Table 5). Apart from this observation, no definite trend was found across all simulations. However, it is observed that change in the SASA value of hydrophobic residues which result in exposure or burial is more often observed in twin-range cut-off simulations. Burial or exposure of hydrophobic residues during the course of PME simulations do not occur frequently. This shows that the dynamics of loop LB in burying the hydrophobic groove is different between the PME simulations and those simulations which used twin-range cut-off.
In summary, the hydrophobic residues in the hydrophobic groove, which were participating in interactions with BH3 peptide ligands, take part in stable interactions among themselves. The  dynamic nature of loop LB connecting the helices H2 and H3 helps to bury up to 250 to 380 Å´2 surface area of bulky hydrophobic residues from exposure to the solvent. Both observations help to explain how the exposed hydrophobic groove remains stable in solvent-exposed environment.

Functional significance of helix H2 unwinding
One of the significant results of apo-and holo-Bcl-X L simulations using twin-range cut-off is the dramatic unwinding of helix H2 which was also observed in the Bcl-X L complex simulations [30]. It is important to note that helix H2 contains the functionally important BH3 domain and such unwinding is not observed in other major helices of Bcl-X L . This observation could have functional significance in two respects and the first one is  These interactions are stable in apo/holo Bcl-X L simulations but not in the Bcl-X L complex simulations [30]. b Residue shown in bold participate in stable interactions with the BH3 peptide ligands in Bcl-X L complex simulations [30]. related to the binding of diverse BH3 peptide ligands to the hydrophobic groove of Bcl-X L . The unfolding of helix H2 can provide flexibility to the hydrophobic groove and hence Bcl-X L can bind to the BH3 regions of different pro-apoptotic proteins with different affinities [5,22]. This is evident from the fact that Bcl-X L interacts with the BH3 peptides of Bak, Bad and Bim proapoptotic proteins and the affinities differ from 0.6 nM to 340 nM [19,21,22]. Structures of Bcl-X L in complex with mutant BH3  SASA was calculated for a given residue using the g_sas tool available in GROMACS version 3.3 [39,40]. b Residues shown in bold participate in stable interactions with at least one of the BH3 peptide ligands in Bcl-X L complex simulations [30]. peptides and BH3 peptidomimetics revealed structural transitions that might take place in the BH3-binding groove [23,49]. Structural plasticity of hydrophobic groove has been demonstrated in other Bcl-2 family structures in complex with BH3 domains of different pro-apoptotic proteins [50]. Plasticity in ligand-binding region has been recognized in other diverse proteins also [51,52,53,54]. Thus the unwinding of H2 can be related to the protein's ability to bind different BH3 domains with differential affinities.
The second aspect of helix H2 unwinding may have implications in the formation of homo/hetero dimers of Bcl-2 members. The diverse Bcl-2 proteins are known to adopt the same helical fold [18,55] and hence the hydrophobic face of the amphipathic BH3 region which is contained in helix H2 will be buried in this fold. In order for the BH3 region of these pro-apoptotic proteins to bind in the hydrophobic groove of their anti-apoptotic partners, the BH3 region has to be exposed. The exposure of BH3 region requires conformational change of helix H2. One suggested change in earlier studies was rotation of helix H2 to expose its hydrophobic surface [19]. However since helix H2 forms part of the protein's hydrophobic groove, a simple rotation with helical structure intact requires breaking of hydrophobic contacts with the rest of the protein and such conformational change may not be energetically feasible. As observed in earlier complex simulations [30] and the present apo-and holo-Bcl-X L simulations, unwinding of helix H2 and later reforming the helix with hydrophobic surface exposed for binding with the hydrophobic groove may be one way that can lead to the protein-protein interactions. Formation of homo-and hetero-dimers has been reported for several Bcl-2 members including Bcl-X L and the interaction occurs when BH3 domain of one protein binds to the hydrophobic groove of another protein [56,57,58]. Several studies have suggested that exposure of BH3 domain is a crucial step for homo-and hetero-dimerization of Bcl-2 proteins. In many studies that investigated dimerization of the pro-apoptotic proteins Bax and Bak, exposure of BH3 domain seems to be critical [35,36,56]. For example, studies on Bak mutants altering its BH3 domain or its hydrophobic groove clearly influenced its ability to form oligomers and this step is required for its pro-apoptotic function [35]. This study demonstrated that exposure of BH3 domain of Bak is an essential early step for its interaction with the hydrophobic groove. The present simulation studies suggest that such exposure of BH3 domains can occur first by unwinding the BH3-containing helix H2 and later reforming the helix by exposing its hydrophobic side to the hydrophobic binding groove. Hence, unwinding of helix H2 can either provide structural plasticity to the hydrophobic groove by allowing diverse BH3 peptides to bind and/or can be part of the process to expose the BH3 domains that can facilitate the formation of homodimers or heterodimers. Either way, this event assumes major functional significance in the apoptotic pathway.
As far as the PME simulations are concerned, unwinding of helix H2 is not observed perhaps due to the fact that the lengths of the simulations are short. A detailed discussion comparing the results of PME and twin-range cut-off simulations is given in a separate section (see below).

Dual role of flexible loop LB
The stability of exposed hydrophobic groove can be explained by two factors: (i) interactions among the hydrophobic residues in the groove and (ii) burial of hydrophobic amino acids by the residues from loop LB. Interactions between exposed hydrophobic residues is akin to hydrophobic collapse when a protein is on the way to its folding process trying to shield its hydrophobic residues from solvent exposure. Hydrophobic residues are also protected from the solvent by loop LB. This loop connects the BH3containing helix H2 with the next helix H3 and it wraps around the hydrophobic groove ( Figure 1). In the apo-Bcl-X L structure, a short stretch of this loop assumes helical conformation while in the Bcl-X L structures in complex with pro-apoptotic BH3 peptides, it essentially adopts a random conformation. It has been shown that the loop LB has stable interactions with the bound BH3 peptides [30]. The present simulation studies clearly demonstrate that this loop plays a major role in shielding the hydrophobic residues in the hydrophobic groove of Bcl-X L . The flexible and dynamic nature of this loop is evident by comparing the apo-Bcl-X L structure with the structures of other Bcl-2 members, structures of Bcl-X L complexes and also the mutant Bcl-X L structures [19,21,31,37,59,60,61]. We have identified 13 structures of Bcl-X L complexes in PDB in which the bound molecules are BH3 peptides or other organic ligands. Superposition of these structures (Figure 7) clearly demonstrates that loop LB displays conformational heterogeneity to accommodate and to have strong interactions with the ligands. In the absence of any bound ligands, this loop helps burying 250 to 380 Å´2 surface area of bulky hydrophobic residues when all eight simulations are considered. The burial and exposure of loop LB residues is not uniform across all simulations. This indicates that loop LB is constantly sampling many conformations to optimally bury as many hydrophobic residues in the hydrophobic groove as possible. Thus loop LB plays a pivotal role in stabilizing the exposed hydrophobic residues in the hydrophobic groove in the absence of any bound ligand. The role of loop LB residues in conferring selectivity and higher affinities to different Bcl-X L ligands has to be explored further.

Twin-range vs PME simulations
Although the twin-range cut-off of 10 to 18 Å´is long enough to include most of long-range interactions, one can always argue that the long-range interactions beyond 18 Å´may still be not negligible. In such a case, the neglected long-range interactions could have influenced the properties and behavior of Bcl-X L protein. Since the development of Ewald-based methods in early 1990s [62,63,64], they have become one of the most popular and widely used methods to calculate the long-range interactions in molecular simulations. Hence in this study, we have also simulated both apoand holo-Bcl-X L in which particle mesh Ewald (PME) scheme was used in calculating the long range interactions. Four simulations (Apo-pme-I, Holo-pme-I, Holo-pme-II and Holo-pme-III) have been carried out as described in the Methods sections and summarized in Table 1. These simulations differed in the initial structure (apo or holo) and the size of the box. In all four simulations, the same equilibration scheme (Scheme I) was used. The systems were simulated for a period of 25 to 55 ns.
When we compared the set of simulations that used twin-rangecut-off with the PME simulations, the most significant difference between the two sets is the helix H2 unwinding in the former while the same helix is very stable in the latter set of simulations. This gives rise to a set of questions. Is such an unwinding due to the neglect of long-range interactions in twin-range cut-off simulations? If it is so, why do most of the other major helices in apo-and holo-Bcl-X L remain stable? Is this the only difference between the twin-range cut-off simulations and PME simulations? To answer the above questions, we first calculated the RMSD profiles of all the simulations with only helices H1, H3, H4 and H5 and the MD trajectories were compared between the two sets of simulations ( Figure 8A and 8B). It is very clear that the RMSD profiles are very similar for the two sets of simulations if helix H2 and the Cterminal H6 are excluded from the analysis. This also accounts for the difference observed in the overall Bcl-X L helical content between twin-range cut-off and PME simulations.
Molecular simulations using PME scheme to calculate longrange interactions have been compared with other methods like reaction field, twin-range cut-off and continuum models [65,66,67,68,69,70,71,72,73,74,75]. While these earlier studies reported less conformational sampling and reduced flexibility in PME simulations, several recent simulation studies have used PME and investigated flexible regions and conformational changes in diverse proteins [76,77,78,79]. Hence, the highly stable helix H2 in apo-and holo-forms in the set of PME simulations in the present study indicates that the length of the simulations is not sufficient to observe the H2 destabilization. The time-period of present simulations, 25 to 55 ns, is perhaps not long enough for the helix H2 to unwind. This is further corroborated by a recent study of Bcl-X L using Amber 99SB force field [80] and PME scheme. This study by Yang and Wang [81] reported only minor backbone changes for the apo-Bcl-X L in water during the 32 ns simulation. Only in the presence of isopropanal cosolvent molecules, large conformational changes were observed in specific regions.
To further investigate the influence of the scheme that is used to calculate the long-range interactions on the helix H2 stability, we carried out six additional simulations ( Table 6) using GROMACS version 4.5.5 [82]. Three simulations used longer twin-range cutoff. In the other three simulations, the truncated helix H6 was extended using the structure of Bcl-X L complex (PDB ID: 1G5J) as a template. Within 10 to 35 ns, unwinding of helix H2 was observed in all the simulations (Figure 9). The choice of simulation conditions can affect the dynamics of the system as was demonstrated recently [83]. Hence, in the set of PME simulations, the dynamic behavior of Bcl-X L showed differences with the twinrange cut-off simulations and we believe that a longer PME simulation would eventually produce a similar destabilization of helix H2.

Conclusions
The anti-apoptotic Bcl-X L protein is an attractive target for anti-cancer drugs. The pronounced hydrophobic groove formed by the helix bundle structure is the ligand-binding region and the Bcl-X L inhibitors are developed keeping the physical and chemical nature of this groove in mind. In the absence of any bound ligand, Figure 8. Comparison of two sets of simulations using RMSD analysis. Comparison of MD trajectories of RMSD values calculated for simulations that used (A) twin-range cut-off and (B) PME to calculate the long-range interactions. Only the helices H1, H3, H4 and H5 were considered for calculating RMSD between each MD simulated structure and the starting Bcl-X L structure. Helix H2 and the C-terminal helix H6 were excluded in this analysis. Compare this figure with Figure 2B and Figure 3B. doi:10.1371/journal.pone.0054397.g008 the present simulation study investigated the behavior of the solvent-exposed hydrophobic groove. Eight independent simulations of Bcl-X L in apo and holo form were carried out. These simulations used either twin-range cut-off or PME for calculating the long-range interactions. All the twin-range cut-off simulations exhibited destabilization of the BH3 domain-containing helix H2. However, the other major helices with the exception of C-terminal helix H6 were stable. At this point, it is speculated that the failure of PME simulations to destabilize the helix H2 may simply be due to the length of simulations. If the PME simulations are extended further, it is possible to detect the conformational changes associated with helix H2. Based on several experimental results, the unwinding of helix H2 can be linked to the plasticity of the hydrophobic groove which enables the Bcl-X L protein to bind to different BH3 ligands with differential affinities. Helix H2 destabilization can also be connected to the formation of homoor hetero-dimers of Bcl-2 proteins. Since helix H2 contains BH3 domain, it has to undergo conformational changes to expose the buried hydrophobic side of BH3 domain. Hence the loss of helical character for H2 seems to have functional significance. Results of simulation studies show that the exposed hydrophobic residues from the groove interact among themselves while in the structures of complexes, most of them were involved in interacting with the BH3 peptide ligands. The solvent accessible surface areas of these residues are significantly buried by the loop LB connecting the helices H2 and H3. This explains how the predominantly hydrophobic groove remains stable when exposed to the solvent. Understanding of the plasticity of hydrophobic groove and the dynamics of loop LB reported in this study can help in the design of inhibitor molecules that will be highly specific to Bcl-X L protein. Helices and side-chains of hydrophobic residues are displayed in ribbon and stick representation respectively. Surface and ribbon representations of helices H2, H3, H4, H5 and loop LD (cyan) along with the hydrophobic residues from these regions (yellow) are shown for (C and E) Apo-pme and (D and F) Holo-pme-I simulations without loop LB (C and D) and with loop LB (E and F). Loop LB surface is represented in purple color in (E) and (F). The Bcl-X L structures shown in this figures were saved at the end of 55 ns production run from Apo-pme and 50 ns production run from Holo-pme-I simulation.

Supporting Information
(JPG) Figure S4 Hydrophobic residues in the hydrophobic cleft: Interactions and accessible surface areas in Holopme-II and Holo-pme-III simulations. Interactions among the hydrophobic residues in the hydrophobic groove are shown for (A) Holo-pme-II and (B) Holo-pme-III simulations. Helices and side-chains of hydrophobic residues are displayed in ribbon and stick representation respectively. Surface and ribbon representations of helices H2, H3, H4, H5 and loop LD (cyan) along with the hydrophobic residues from these regions (yellow) are shown for (C and E) Holo-pme-II and (D and F) Holo-pme-III simulations without loop LB (C and D) and with loop LB (E and F). Loop LB surface is represented in purple color in (E) and (F). The Bcl-X L structures shown in this figures were saved at the end of 25 ns production runs from Holo-pme-II and Holo-pme-III simulations. (JPG) Apo-H6-extended-II Simulation was started with different initial velocities using Apo-H6-extended-I 10 ns Apo-H6-extended-III Positional restraints were applied on heavy atoms of helix H6 c 10 ns a In all listed simulations, Scheme I equilibration protocol was followed. All other simulation parameters are provided in the Methods section. b Bcl-X L structure with PDB ID 1G5J was used as a template to model the extended helix H6. Modeller ver 9.11 [85]