A study of the structural properties and thermal stability of human Bcl-2 by molecular dynamics simulations

The anti-apoptotic B-cell lymphoma 2 (Bcl-2) protein interacts with several proteins that regulate the apoptotic properties of cells. In this research, we conduct several all-atom molecular dynamics (MD) simulations under high-temperature unfolding conditions, from 400 to 800 K, for 25 ns. These simulations were performed using a model of an engineered Bcl-2 human protein (Bcl-2-Δ22Σ3), which lacks 22 C-terminal residues of the transmembrane domain. The aim of this study is to gain insight into the structural behavior of Bcl-2-Δ22Σ3 by mapping the conformational movements involved in Bcl-2 stability and its biological function. To build a Bcl-2-Δ22Σ3 three-dimensional model, the protein core was built by homology modeling and the flexible loop domain (FLD, residues 33-91) by ab initio methods. Further, the entire protein model was refined by MD simulations. Afterwards, the production MD simulations showed that the FLD at 400 and 500 K has several conformations reaching into the protein core, whereas at 600 K some of the alpha-helices were lost. At 800 K, the Bcl-2 core is destabilized suggesting a possible mechanism for protein unfolding, where the alpha helices 1 and 6 were the most stable, and a reduction in the number of hydrogen bonds initially occurs. In conclusion, the structural changes and the internal protein interactions suggest that the core and the FLD are crucial components of Bcl-2 in its function of regulate ng access to the recognition sites of kinases and caspases.


Introduction
Apoptosis is a natural process required for the removal of redundant cells during development, potentially dangerous cells and cells in senescence (Kerr, Wyllie, & Currie, 1972). Thus, the cell death dysregulation has been implicated in a variety of human diseases such as cancer (Filhol & Cochet, 2011), autoimmunity disorders (Otsuki et al., 2011), and neurodegenerative disorders (Ola, Nawaz, & Ahsan, 2011). This process is regulated by several proteins that belong to the Bcl-2 family, which have been extensively studied. Members of this family are grouped according to their homology and participation in the apoptotic mitochondrial pathway (proand anti-apoptotic proteins) (Chipuk & Green, 2008). The anti-apoptotic Bcl-2 proteins (e.g. Bcl-2, Bcl-xL and Bcl-W) possess up to four Bcl-2 homology (BH) sequence domains (BH1-BH4), whereas pro-apoptotic proteins (e.g. Bax, Bak, and Bok) have either three BH domains (BH1-BH3) or only the BH3 domain (e.g. Bid, Bin, and Bik) . It is known that BH domains are involved in protein interactions among Bcl-2 family members. However, these interactions at the 3D structural level are scarcely studied, and consequently, there is not enough structural information about the proteins' behavior during either anti-or pro-apoptotic events (Lin et al., 2004;Oltvai, Milliman, & Korsmeyer, 1993). There is experimental evidence showing that this protein family undergoes structural and/or activity changes during phosphorylation (Fang et al., 1998;Yamamoto, Ichijo, & Korsmeyer, 1999) and proteolytic cleavage (Cheng et al., 1997;Clem et al., 1998). Particularly, the Bcl-2 protein was identified originally at the breakpoint of a t(14:18) translocation in a lymphocytic leukemia cell line (Tsujimoto, Finger, Yunis, Nowell, & Croce, 1984). In addition, the Bcl-2 protein has been found to be over-expressed in many cancer cells such as B cell derivative lymphomas (Son et al., 2010) and colorectal adenocarcinomas (Lewis, Thrumurthy, Pritchard, Armstrong, & Attwood, 2011). Additionally, the Bcl-2 protein has been involved in the chemotherapeutic resistance of many cancers (Deng, Kornblau, Ruvolo, & May, 2001;Wesarg et al., 2007). For that reason, some work-groups have focused on designing new compounds to inhibit Bcl-2 activity .
Computational studies are of great interest for describing the folding/unfolding properties of proteins and for explaining the interaction properties with other proteins or small ligands, which can be used in rational drug design. Additionally, active/inactive structural conformations can be correlated with experimental data (Budi, Legge, Treutlein, & Yarovsky, 2004). Although there are several experimental assays that can achieve protein unfolding (Bennion & Daggett, 2004;Carra & Privalov, 1996;Chan, Bromberg, & Dill, 1995), molecular dynamics (MD) simulations allow the exploration of the protein model while considering various biological properties (Wang et al., 2010). Extreme thermal conditions in a MD simulation can accelerate the protein conformational movements, giving insight into the unfolding properties and an atomistic level description of these events (Daggett & Levitt, 1994;Gu, Wang, Zhu, Shi, & Liu, 2003). Recently, experimental results have suggested several hypotheses about the behavior of the Bcl-2 protein (Blagosklonny, 2001;Chang et al., 1997;Deng et al., 2006;Haldar et al., 1998;Huang et al., 1998;Kazi et al., 2011;Kim et al., 2001;Kirsch et al., 1999;Ku et al., 2011;Lin et al., 2004;Petros et al., 2004;Reed et al., 1996;Ruvolo et al., 2001;Srivastava et al., 1999). However, to our knowledge, this protein has not been studied using MD simulations at different high-temperature conditions to describe its structural properties. A similar effort was reported by Raghav et al. where MD simulations were used to study the flexibility and stability of the FLD, as well as its binding properties to other related proteins (Raghav, Verma, & Gangenahalli, 2012). Other MD simulations of Bcl-2 family proteins such as Bax, Bcl-xL, Mcl-1, and BNip3 have been performed using diverse force fields (Rosas-Trigueros, Ilizaliturri-Flores, Benítez-Cardoza, Correa-Basurto, & Zamorano-Carrillo, 2012). A key functional domain in the Bcl-2 family proteins is represented by BH3 because it binds with less selectivity, compared with other domains, and is critical for the survival of the cell (Smits, Czabotar, Hinds, & Day, 2008). MD simulations on BH3 domain have proposed that the stability of its helical character depends on the length of each helix and their interaction with water. Interestingly, it has been proposed that there might be a correlation between the stability of the BH3 domain and its affinity for Bcl-xL (Bharatham, Chi, & Yoon, 2011). Complexes of Bcl-2 family proteins or regions of these such as Bcl-xL/Bim or Bcl-xL/Bak(BH3) have been also simulated. It has been reported that the stability of the secondary structure might be also been associated to the stability of the interactions (Lama & Sankararamakrishnan, 2011;Pinto, Orzaez, Delgado-Soler, Perez, & Rubio-Martinez, 2011). The dynamics of the FLD in Bcl-2 that contains a sequence of residues described as intrinsically disordered region (IDR) has been characterized in terms of stability and flexibility by MD simulations (Hinds et al., 2007). This functional combination of folding and binding can be provided by IDPs, in particular to the anti-and pro-apoptotic proteins, offering a trait of plasticity to the interactions that could be finely tuned by a multitude of conformations that dynamically change over time and population (Koonin & Aravind, 2002;Moroy, Martin, Dejaegere, & Stote, 2009;Romero et al., 2006). In this work, we carry out homology modeling and an ab initio procedure to obtain the Bcl-2 structure, including its FLD, followed by refinement under MD simulations. Later, the whole system was heated at different temperatures and several structural changes were assessed. The 3D model used in these MD simulations is inspired by the functional Bcl-2 protein expressed and described by Kim et al. (2001), and also three N-terminal residues were include, name here Bcl-2-Δ22Σ3.

Theoretical procedure
To obtain the initial coordinates of these simulations, we constructed a 3D model of the Bcl-2 sequence lacking 22 C-terminal residues . Additionally, three residues in the N-terminal region (Gly1, Ser2, and His3) were added to compare this model with a recombinant protein produced in our laboratory group (unpublished data), since these three N-terminal residues were added to the construct during the cloning process. Bcl-2-Δ22Σ3 sequence was submitted to the I-TASSER server (Zhang, 2008) which used homology modeling for the Bcl-2 protein core and ab initio for FLD, the NMR structure of Bcl-2 (PDB code: 1GJH) was used as template. The strategy of I-TASSER server includes a modeling of fragments that evade being modeled by threading, i.e. unaligned regions that difficult the reassembling of the full-length model. These regions, mainly loops, are built by ab initio modeling, where the peptide is considered as a lattice chain connecting C(alpha) atoms, with attached C(beta) atoms and side-chain centers of mass. Also, the interactions are governed by a force field that includes short-range and long-range knowledge-based potentials obtained from a statistical analysis of the properties of protein structures (Zhang, Kolinski, & Skolnick, 2003). The Ramachandran plot of this model indicated that none of the non-Gly residuepoints fall in disallowed regions (data not shown). In addition, a structural analysis of the model was achieved by the validation programs VERIFY3D and Errat servers, confirming that our model is suitable to be used to provide the initial coordinates for the rest of the MD simulations. Then, hydrogen atoms were added using the software "psfgen" from the visual molecular dynamics (VMD) program (Humphrey, Dalke, & Schulten, 1996). Periodic boundary conditions and electrostatic interactions were calculated by the particle-mesh Ewald method (Batcho, Case, & Schlick, 2001). All MD simulations were performed with the NAMD program 2.8 (Phillips et al., 2005) employing the CHARMM 27 force field (MacKerell et al., 1998). Non-bonded interactions were calculated for a cut-off of 10 Å with a switching function to 8 Å. The list of non-bonded interactions included atoms at a distance less than 11.5 Å. The SHAKE method (Ryckaert, Ciccotti, & Berendsen, 1977) was used with an integration step of 2 fs, while keeping all bonds of the hydrogen atoms rigid. The Bcl-2-Δ22Σ3 3D model was submitted for energy minimization in vacuo. The protein was then solvated with 12,719 water molecules, and one sodium ion was added to equilibrate the electric charges making a rectangular box system (TIP3P; $80 Â 110 Â 78 Å 3 ). This entire system was minimized and treated under MD simulations, maintaining the fixed protein backbone in both procedures. The minimization process consisted of three stages of 1500 steps, each stage using the conjugate gradient algorithm. After these preparatory steps, the system was heated from 0 K to the desired temperature (310, 400, 500, 600, and 800 K) in the first 60 ps of the 100 ps that constituted the equilibration period without restrictions. Then, MD simulations were performed using an NVT ensemble during 25 ns, and the trajectory was stored every 1 ps. The results were analyzed using the Carma program (Glykos, 2006) to obtain the root mean square deviation (RMSD), radius of gyration (Rg), protein solvent accessible surface area (SASA), and root mean square fluctuation (RMSF). The VMD program was used to analyze all the secondary and tertiary structures.

Bcl-2 conformational changes
The Bcl-2 protein has been found to be over-expressed in many cancer cells and has also been shown to be related to the resistance of these cells to anticancer drugs (Masood, Azmi, & Mohammad, 2011). For this reason, this protein has been proposed as a target for drug design. Some reports have described the Bcl-2 FLD as having residues susceptible to phosphorylation, as well as cleavage sites for caspase-3 that could regulate the activity of the protein (Blagosklonny, 2001;Kirsch et al., 1999). In addition, the FLD has not been studied at the experimental level because it has a very high B-factor due to its great mobility. Therefore, this FLD could be considered as an IDR (Raghav et al., 2012). To obtain insights related to Bcl-2 function and its structural properties, MD simulations that reproduce and explain the experimental results are required. Regarding the Bcl-2 protein, the presence of the transmembrane domain makes the expression of the whole protein by heterologous systems difficult. Additionally, the structural studies have been hampered by the poor solubility of the protein. Consequently, some mutations and structural modifications were required to obtain the recombinant Bcl-2 for in vitro studies . On these grounds, we used a reported incomplete sequence to obtain a 3D model of Bcl-2 ( Figure 1(A)). This sequence Bcl-2-Δ22Σ3 was submitted to the well-known homology modeling server I-TASSER, yielding a 3D model of the whole protein combining several folding protein principles (comparative modeling, ab initio folding and threading methods) to obtain the complete 3D model (Roy, Kucukural, & Zhang, 2010). The Ramachandran plot of this model indicated that none of the non-Gly residuepoints fall in the disallowed regions; and RMSD values, regarding alpha carbons between generated model and a reported Bcl-2 structure (PDB:1GJH) was 1.3 Å. The Bcl-2-Δ22Σ3 3D model that was selected presents a tail (residues 206-220), a core (residues 10-32 and 92-205) and a FLD (residues 33-91) showing an amphipathic alpha-helix (residues 68-74) named αL (Figure 1(B)). These results are in agreement with reported data (PDB code: 1GJH) used as template (Figure 1(B)). In addition, the contact map (cut-off 10 Å) and matrix distance of the structure Bcl-2-Δ22Σ3 was built to examine its interactions between its domains (Figure 1(C)). All these analyses provide confidence that the 3D model presents a high-quality structure that is suitable for performing MD simulations with the purpose of refining the structures and submitting the modeled protein to heat conditions.

MD simulation at 310 K
We explored the conformational space of Bcl-2-Δ22Σ3 around a local energetic minimum by a MD simulation at 310 K. As it is shown in Figure 2(A), the FLD has high mobility in comparison with the rest of the protein, which explains the difficulties in structurally characterizing this domain. It was observed that the eight alpha-helices are maintained during all MD simulation times (25 ns), and a helix was formed at the amino terminus (residues 1-4), which was stable only during the first 10 ns (Figure 2(B)).
To study the global behavior of the protein during simulations, we assessed several structural parameters.
The RMSD values converged approximately in 2 ns (Figure 3(A)), which are in good agreement with literature reports [41], in like manner the Rg values reached 22.5 (±1.3) Å at 6.6 ns (Figure 3(B)). Despite the great structural movements of the FLD, the RMSD, and Rg values of the whole protein converged without significant changes. Additionally, the initial number of hydrogen bonds (less than 40) increased to 50 at 310 K in the MD simulations (Figure 3(E)). These results indicate that the Bcl-2 structural changes were not severely affected by the mobility of the FLD (Figure 2(A)). In addition, the SASA value remains near 15,000 Å 2 (Figure 3(C)). Relevant movements of the secondary-structural components or domains are illustrated by RMSF plots. In particular, residues 143-158, corresponding to BH1, presented the lowest RMSF values (Figure 3(D)) indicating that this secondary structure is very stable and might support several stress conditions before unfolding. In contrast, the most mobile regions (data not shown) were those with non-regular secondary structure (loops) highlighting the FLD, as expected. It is known that many proteins undergo conformational changes during the performance of their functions, wherefore flexibility is required to achieve biological activity (Rosales-Hernández, Mendieta-Wejebe, Trujillo-Ferrara, & Correa-Basurto, 2010). To give insight into physiological structural changes of Bcl-2, MD simulations are useful for describing these conformational movements at the atomic level. Because the Bcl-2 FLD has not been solved experimentally due to its large structural movements, the analysis of this region is necessary to clear up some structural details. The MD simulations at 310 K showed that some neighboring residues within the FLD interact with each other to form a short alpha-helix (αL) that could be important for making protein interactions. Ser73, located in this short alpha-helix, is susceptible to phosphorylation. Interestingly, Ser73 phosphorylation inhibits the anti-apoptotic function of Bcl-2, preventing hetorodimerization with Bax and thereby inducing apoptosis (Yamamoto et al., 1999).

MD simulations at 400-800 K
As it was shown in the previous section, Bcl-2 undergoes minor structural changes at 310 K (data not shown). However, when high-temperature conditions  were applied to the system, the protein displayed larger conformational changes. At 400 K, it was evident that the protein maintained its native secondary structure, and the FLD was the most dynamic region moving towards the protein core (Figure 4, first row). The FLD interaction with the protein core was maintained from 2.5 to 7.3 ns. Afterwards, the Bcl-2 FLD separated from the protein core at 9.4 ns, but at 22 ns the FLD again reached completely into the Bcl-2 core and maintains this position up to 25 ns. This conformational change suggests a regulatory process by the FLD to block interacting proteins from accessing the BH3 domain. Also, a possible stabilizing interaction was detected, the salt bridge Arg71-Asp105 ( Figure 1S, see supplementary material). Other relevant interactions suggested by the simulations at 400 K from 2.5 to 7.4 ns were Lys20-Ala83, Pro78-Arg101, and Ala79-Arg101. A similar strategy using MD simulations at high-temperature conditions has been used to explore the switch between active/inactive conformations for other proteins like superoxide dismutase (Camperchioli et al., 2011).
At 500 K, shown in Figures 4 and 5(B), the destabilization of the secondary structure was more evident after 10 ns, but the internal alpha-helices were largely maintained. Similarly, at this temperature, as it was shown at 400 K, the FLD reaches the protein core (from 2.8 ns) maintaining this conformation for the rest of the MD simulations. This result means that high-temperature conditions represent a physical stress that can promote structural changes of the FLD, so that this region can reach the protein core, influencing Bcl-2 recognition and interactions with other proteins. However, during the last ns of the MD simulations, the Bcl-2 core was destabilized and the internal secondary structural elements moved away from each other (Figures 4 and 5(B)). These structural modifications are in good agreement with the changes in the number of hydrogen bonds (Figure 3(E)); the highest number of hydrogen bonds is 41 at 20 ns and decreases at the final snapshot to 37. As it is shown in Figure 5(C), the majority of alpha-helices of Bcl-2 are lost in the last 5 ns of the simulation at 600 K, although alpha-helices 2 (corresponding to BH4) and 7 appeared to be more stable. In addition, at 600 K the global structure of Bcl-2 alternates between extensive and compact forms, mainly involving the FLD. Particularly, at 12.8 ns, the Bcl-2 core lost internal interactions (Figure 4). In addition, the number of hydrogen bonds is reduced (approximately 25, 12, and 17 ns) (Figure 3(E)) as a consequence of this unfolding.
Simulations at 800 K (the most extreme temperature tested in this work) included the most severe structural changes. Therefore, rather than showing functional information, the ensemble obtained at this temperature provides information about the stability and the unfolding pathway in the context of atomic motions. Thus, before 5.4 ns, the global structure of the Bcl-2 can still be visualized as two parts: the FLD and the rest of the protein. This separation later terminates, including a loss of compactness of the core. During the rest of the simulation at this temperature, the Bcl-2 core is highly disordered. The secondary structure at 800 K is almost absent except for alpha-helix 1 and, to a lesser extent, alpha-helix 6 (Figures 4 and 5).

BCL-2-Δ22Σ3 RMSD values
The RMSD parameters permits quantification of the structural changes from initial coordinates to subsequent structures sampled during the simulation. At 310, 400, and 500 K, the RMSD values of Bcl-2 converged at approximately 2, 6, and 4 ns, respectively. However, at 600 and 800 K, the rapid increment of the RMSD is terminated at 4 ns; however, fluctuations (of approximately 5 Å) continue to occur (Figure 3(A)). According to the level of thermal stress, diverse conformations of Bcl-2 are obtained. In particular, simulations at temperatures less than or equal to 500 K provide possible structures related to active or inactive forms, whereas at temperatures above 500 K, an unfolding process could be proposed (Chen, Li, & Ma, 2009;Rosas-Trigueros, Correa-Basurto, Benítez-Cardoza, & Zamorano-Carrillo, 2011). Concerning functionality, snapshots show that the conformational movements led to the interaction between the Bcl-2 FLD and the protein core (Figure 4), which might be related to a regulation of Bcl-2 phosphorylation, caspase cleavage, and the interaction of Bcl-2 with pro-apoptotic proteins (Blagosklonny, 2001;Deng et al., 2006;Haldar et al., 1998;Huang et al., 1998;Kazi et al., 2011;Kirsch et al., 1999;Lin et al., 2004;Reed et al., 1996;Ruvolo et al., 2001;Srivastava et al., 1999). Regarding unfolding, this work proposes alpha-helices 1 and 6 as the most stable secondary structural elements, notwithstanding RMSD fluctuations.

BCL-2-Δ22Σ3 hydrogen bonds
The hydrogen bonds of Bcl-2 are involved in internal protein interactions such as residue-residue interactions mediated by the side-chain or backbone, or inter secondary structures. The number of hydrogen bonds decreases when increasing the temperature of the simulation (shown in Figure 3(E)). For instance, at 310 K, we observed an average of 45 hydrogen bonds, while at 800 K less than 20 hydrogen bonds remained. Snapshots at these temperatures allow the observation of the lost internal contacts within the Bcl-2 protein, as well as the extension of the protein (Figure 4). Experimental research has shown a decrease in internal hydrogen bonds occurs during thermal denaturation processes (Mancinelli, Caraglia, Budillon, Abbruzzese, & Bismuto, 2006).

BCL-2-Δ22Σ3 Rg values
Additional information about Bcl-2 mobility/flexibility was obtained by analyzing the time course of the Rg profile along the MD simulations at different temperatures. This parameter quantifies the proximity of each atom to the center of mass of the protein. At 310 K the Rg values along the MD simulation are between 22 and 24 Å. This distance can be explained by the presence of the extended FLD protruding away from the protein core. In turn, the Rg values are reduced when the temperature increases: no more than 20 Å at 400 and 500 K, and values from 18 to 22 Å at 600 K following to the movement of the FLD towards the protein core. On the other hand, at 800 K a Rg value of 20 Å (2 ns) is the starting point for a continuous increase until 24 Å, approximately. Nevertheless, three prominent peaks at 8, 16, and 19 ns suggest drastic increments in the unfolding of the Bcl-2 followed by a subsequent compaction (Figure 3(B)).

BCL-2-Δ22Σ3 RMSF values
To detect the most mobile domains in the Bcl-2 protein, the RMSF of each residue was calculated. Specifically, we evaluated the backbone-residue movements at different temperatures, which are shown in the RMSF plot.
The FLD and C-terminus had marked fluctuations at 310, 400, and 500 K. Even at 600 K, the FLD represents, together with its vicinity (residues 22-122), the largest contribution. This high mobility might explain the difficulty in solving the FLD structure (PDB code: 1GJH).
As expected, the structures generated at 800 K showed higher RMSF values. In particular, a region located between residues 156 and 172 (helix 6) underwent severe fluctuations (Figure 3(D)).

Identities and differences between Bcl-2 and Bcl-xL
It is known that Bcl-2 and Bcl-xL show the same antiapoptotic function. Consequently, their expression has been detected in common tumors, and common inhibitors can bind both proteins (Arbab et al., 2012;Oakes et al., 2012). The similar structure and properties of Bcl2 and Bcl-xL allow to build the Bcl-2 core by including a short loop of Bcl-xL instead of the FLD (PDB: 1GJH and 1G5 M). However, there are notable differences that lead to the recognition by distinct peptides (Petros et al., 2001). Thus, some residue protein sequences consequently could explain their differential affinity for ligands. The study of the Bcl-2 structural behavior under MD simulations, taking into account the FLD, suggests the mechanistic of its functional domains. Hypothesizing, whether, some of these mechanisms are shared or not with other Bcl-2 family members is a matter for further experimental and computational studies.

Functional implications of Bcl-2 movements
The structural movements of the Bcl-2 protein depict some insights that can be associated with the Bcl-2 protein function. To this end, the structure-functionmovement relationship in this protein is explored by MD simulation under thermal stress of a monomer. Interestingly, the Bcl-2 protein has a BH3 domain that is not able to reach the hydrophobic groove of another Bcl-2 as occurs in pro-apoptotic proteins (Borner et al., 1994;Dlugosz et al., 2006). However, there is evidence that in the absence of the FLD, the BH3 domain is exposed and prone to interact with the BH3 of another Bcl-2. Consequently, Bcl-2 can act as an anti-apoptotic protein (Chang et al., 1997;Cheng et al., 1997). The FLD also has very important regions that are implicated in the Bcl-2 inactivation following Ser73 phosphorylation (Yamamoto et al., 1999). The results of the MD simulations at 400 K show that the Bcl-2 FLD reaches into the protein core and interacts mainly with the BH3 domain. This might have important implications in understanding the swapping between active/inactive conformations of Bcl-2. The FLD does not stay blocking the BH3 domain over the whole time course of MD simulations; instead, it acts as a lid that opens and closes towards the BH3 domain, i.e. switching between possible active and inactive forms of Bcl-2. This switch might be triggered by the phosphorylation of Ser73 (located in the αL helix), regulating the movement of the FLD  to interact with the Bcl-2 core (Deng et al., 2006). Thus, conformational movements induced by high temperatures suggest that to avoid interactions of Bcl-2 with pro-apoptotic proteins, the FLD hides the core of the protein. However, there is no experimental evidence of this observation (Boukharta et al., 2011). Another reported functional role of the FLD is its participation in regulating cell survival by interaction with p53, a key cell cycle arrest protein activated in response to DNA damage (Deng et al., 2006). The disposition of the FLD domain could favor a suitable interaction between Bcl-2 and p53.
MD simulations also provided certain insights about the C-terminus, which plays a crucial role as a metabolic sensor (Azad et al., 2006). During the MD simulations at 400 K, the C-terminus reaches the core and the BH3 domain, which could have influence on the inactive properties of the Bcl-2 protein. It is necessary to note that the C-terminus is anchored to the mitochondrial membrane and is able to interact with the hydrophobic Bcl-2 core when it is not inserted into the membrane (Szilák, Moitra, Krylov, & Vinson, 1997). In addition, the anti-apoptotic function of Bcl-2 is independent of anchoring (Moroy et al., 2009). Other experimental evidence shows that during the apoptosis, Cys158 moves from a hydrophilic to a hydrophobic environment (Kolmar, 2011). This movement might be correlated with the high mobility shown in this amino acid, according to the RMSFs at 500 and 800 K, explaining the possible role of Cys158 in triggering apoptosis  and its association with the mitochondrial membrane (Szilák et al., 1997).

Stability of the native conformation of Bcl-2 at different temperatures and folding/unfolding mechanisms
Protein stability is still an open problem, and computational models provide information to address it, emphasizing their application in drug design (Arbab et al., 2012;Baldwin, 1975;Onuchic & Wolynes, 2004). In the case of Bcl-2, part of the secondary structure was preserved even at 600 K, traces of alpha-helices were detected ( Figure 5(C)). Meanwhile, at this temperature, the tertiary structure has undergone several changes as shown by geometric parameters (Figure 3). However, with up to 800 K of temperature applied to the system, Bcl-2 exhibited a larger extension. The high thermal stability might be explained by the internal non-bonded interactions found in Bcl-2 as occurs in Bax protein, a representative pro-apoptotic member of the Bcl-2 family (Ko et al., 2011).
Several structural details that might be implicated in Bcl-2 functionality were observed at 500 K. At this temperature, Bcl-2 shows greater conformational changes, but the secondary structure is mostly maintained and Bcl-2 resists unfolding ( Figures 4 and 5(B)). At 600 K, Bcl-2 shows larger solvated surfaces due to the extended conformations in which the protein chain is extensively hydrated, and the individual residues do not interact with each other as in the native form (Figure 3 (C)). Based on the knowledge acquired from experimental and computational studies, several models have been proposed for protein unfolding (Onuchic & Wolynes, 2004). Interestingly, it has been demonstrated that, for some proteins, the hydrophobic core is readily structured before the formation of secondary structures (Chan et al., 1995). In contrast, the framework model states that formation of the tertiary contacts is preceded by the completion of the secondary structures (Baldwin, 1975). These data suggest that in the case of Bcl-2, the framework model might be more suitable to describe the folding/unfolding reaction of Bcl-2.

Concluding remarks
A Bcl-2 3D model including the FLD was obtained by homology modeling and ab initio methods. Furthermore, the conformational space and stability of the Bcl-2 protein were explored by MD simulations at high-temperature conditions, and several structural parameters were analyzed. Bcl-2 shows an overall high stability because it maintains its globular form despite up to 600 K of heating. Interestingly, these results indicate that, at 400 and 500 K, the FLD comes close to the protein core interacting with the BH3 domains, avoiding exposure to possible interactions with other members of the Bcl-2 family. The hydrophobic face of helix αL is also oriented toward the Bcl-2 core. Finally, the secondary structures at 600 and 800 K show that alpha-helix 1 and 6 are the most stable elements. These findings improve the understanding of the stability of Bcl-2 and the switching between the active/inactive forms of this protein to design drugs that regulate the cell death.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.10.1080/07391102.2013.833858.