Molecular Determinants for the Activation/Inhibition of Bak Protein by BH3 Peptides

Apoptosis is a key procedure for all cells. Understanding this process and its regulation has been a subject of study in the last decades. Bcl-2 family of proteins are involved in the regulation of the apoptosis through the formation of homodimers or heterodimers between anti-apoptotic and pro-apoptotic proteins. Deregulation of pro-apoptotic proteins contributes to the progression of many tumour processes. Understanding how these proteins are activated is key to find new anti-cancer treatments. As no drug capable of activating Bak has been discovered yet, studying the structural and energetic insights of the binding of the known Bak activators, BH3 peptides, is essential to design new small molecules that resemble their binding to Bak. Recently, a BH3 Bim analogue has been discovered which inactivates Bak instead of activating it. Therefore, the present work is aimed to identify how BH3 peptides activate or inactivate Bak and determine any difference between them. Determination of the structural differences between the complexes with the activator and the inhibitor has been carried out by means of the study of the fluctuations of the corresponding Principal Components. Moreover, to calculate the binding free energy of the different complexes and to determine which residues of the peptide have the largest contribution to complex formation, MMPB/GBSA approaches are used. Results obtained in this work show differences in complexes with the activator and the inhibitor in structural and energetic terms, which can be used in the design of new molecules that could activate or inactivate pro-apoptotic Bak.


Introduction
Apoptosis is a form of programmed cell death indispensable for tissue homeostasis, embryonic development and immunity. 1 In this process, a cell undergoes a series of morphological transformations that lead to phagocytes recognition with its subsequent engulfment. 2 Deregulation of apoptosis has a major impact in disease. Thus, an excessive response to apoptosis signals can lead to an increased ischemic condition or drive neurodegeneration, whereas defective apoptosis has a major role in tumour development and autoimmune diseases, among others. 3 Apoptosis is orchestrated by diverse cysteine proteases of the caspase family 4 and is initiated through two different signalling pathways: the extrinsic and intrinsic pathways. 1 In the former, the process is mediated by cell surface death receptors that recognize specific signals from other cells and activate the caspase cascade, whereas in the latter, mitochondrial outer membrane permeabilization (MOMP) is produced as a consequence of cell stress or damage. 1,5 MOMP facilitates the release of pro-apoptotic factors, such as cytochrome c, into the cytosol, triggering the assembly of the caspaseactivating complex to initiate the activation of the executioner caspases 3, 6 and 7 for dismantling of the cell. 1,5 This is a highly regulated process, since once MOMP is produced even if caspase activity is blocked, cells experience a bioenergetic collapse. 6,7 Members of the B-cell lymphoma-2 (Bcl-2) family of proteins are primary responsible for the control of the apoptotic intrinsic pathway, through the regulation of mitochondrial outer membrane permeability. 8,9 These proteins share one or more specific conserved regions known as Bcl-2 homology (BH) domains, recognised to be crucial for function, as deletion of these domains via molecular cloning affects survival/apoptosis rates. Taking into account the number of BH domains they exhibit and their function in the apoptotic process (either pro-apoptotic or anti-apoptotic), these proteins can be classified into three groups. First, the anti-apoptotic or pro-survival Bcl-2 proteins (Bcl-2, Bcl-xL, Bclw, Mcl-1 and Bfl-1/A-1) that exhibit up to four BH-domains (BH1-BH4) and their role is to inhibit the activity of pro-apoptotic Bcl-2 proteins by binding them. Second, the pro-apoptotic effector Bcl-2 proteins (Bak and Bax), exhibiting also up to four BH domains (BH1-BH4), being the BH3-domain necessarily present for mitochondrial apoptosis signalling. These proteins are involved in mitochondrial outer membranes pore formation through homo-oligomerization. Third, the also proapoptotic BH3-only proteins (Bim, Bid, Puma, Bad, Bik, Bmf, Hrk, Noxa), sharing little sequence homology apart from the BH3-domain. This group can be further divided into two subgroups: the direct activators (Bim, Bid and Puma) that directly activate pro-apoptotic effectors, and the sensitizers (all BH3-only proteins), that indirectly activate Bak and Bax by binding to anti-apoptotic proteins and liberating the pro-apoptotic effectors. Some members of the family have also transmembrane domains at their C-terminus, which primarily function to localize them into the mitochondrial outer membrane. 10 The interaction network between members of the Bcl-2 family is dominated by proteinprotein interactions. Indeed, the elucidation of the molecular structure of diverse members and different BH3-domain-Bcl2 complexes has been key in recent years to understand MOMP regulation and provide novel opportunities for therapeutic intervention. 11,12 All multi-BH domain members exhibit a conserved all-α globular structure that consists of nine α-helices (α1-α9) packed around helix α5. The fold exhibits a core domain defined by helices α2 -α5 and a latch domain, defined by helices α6 -α8 (Figure 1) that is involved in the activation process. In the core domain, helices are packed in such a way that define a hydrophobic surface groove that is occupied by the C-terminus of the protein (α9) when it is in its cytosolic form, but also is the site where the BH3-domain of other family members bind to. 13 In contrast, most BH3-only proteins are intrinsically disordered, binding through the BH3-domain into the hydrophobic groove of multi-BH domain members. Bid is the exception to this rule, since it exhibits a multi-BH domain fold and requires cleavage to expose its BH3-domain, generating the active form (truncated form or tBid). 14 Bcl-2 proteins interact with each other generating a complicated interaction network that regulates the apoptotic machinery. Despite knowledge accumulated in the last decades, our understanding of the mechanism of this regulation is not yet complete. 15,16 Thus, it is well established that the relative concentrations of pro and anti-apoptotic proteins are a necessary condition that determines cell commitment to undergo apoptosis. Thus, under normal conditions in healthy cells, pro-survival Bcl-2 proteins sequester pro-apoptotic effectors as well as BH3-only proteins, preventing apoptosis. However, upon a cytotoxic stimulus, the overwhelming number of BH3-only proteins produced activate pro-apoptotic effectors either by direct binding or through the liberation of restrained pro-apoptotic effectors by binding to pro-survival Bcl-2 proteins with the subsequent MOMP induction to result in apoptosis. 9,10 However, the variable affinities exhibited by the diverse members of family for each other and their modulation when proteins are embedded in a membrane are also relevant. 10 In healthy cells, the pro-apoptotic effectors Bak and Bax are shuttled between the mitochondria outer membrane and the cytosol at different rates. 17,18 This differential behaviour makes the former to be primarily attached to the mitochondria outer membrane through its C-terminus (α9), whereas the latter is primarily found in the cytosol, adopting a globular form with its transmembrane domain bound to its hydrophobic groove. 17 Upon cytotoxic stress, Bak and Bax are accumulated and activated at the mitochondrial outer membrane with their subsequent homo-oligomerization that produces MOMP, a crucial event of intrinsic apoptosis. Accumulation of Bax in the mitochondrial outer membrane is executed by an overwhelming number of BH3-only proteins produced after a cytotoxic stress, some activating Bax and others, inhibiting pro-survival Bcl-2 proteins.
Crystallographic studies show that Bax exhibits two binding sites: the canonical hydrophobic groove defined between helices α2 -α5 and a rear activation site, found between helices α1 and α6 at the opposite surface of the hydrophobic groove. 19 BH3-only proteins interact with Bax at the rear activation site, promoting the displacement of the loop between helices α1 and α2, which leads to the displacement of its transmembrane domain from the hydrophobic groove, facilitating membrane insertion. 20,21 On the other hand, as Bak is constitutively inserted in the mitochondrial outer membrane (MOM) through its α9 helix, the canonical hydrophobic groove was postulated as the only activation site. 20,22,23 However, it has recently been identified another activation site in the loop between α1 and α2 triggered by antibodies. 24 Despite this apparent differential behaviour, there is enough experimental evidence to support that they present conserved activation mechanisms. 25,26 Thus, the activation mechanism of these pro-apoptotic proteins once they are inserted in the MOM consists in the transient binding of BH3-only proteins like Bim or tBid producing the dissociation of α1 from the α6 helix, and allowing the separation of the latch domain (helices α6-α8) from the core domain (helices α2-α5). At this time, the affinity of the BH3-only protein decreases and Bak/Bax exhibits fully exposed its BH3-domain to bind another Bak/Bax protein to produce symmetric dimers. These dimers associate via α6:α6 interface and further dimer-by-dimer homo-oligomers and form macropores on the membrane triggering apoptosis. 27 The aim of this work is to analyse the structure of pro-apoptotic Bak and provide the molecular determinants from the analysis of the complex with diverse BH3 peptides involved in the binding and activation of Bak, such as Bid (EDIIRNIARHLAQVGDSMDRSI) or Bim (DMRPEIWIAQELRRIGDEFNAYYARR). The molecular determinants of a recently disclosed Bim analogue are also analysed: Bim(h3Pc) (DMRPEIWIAQELRRXGDEFNAYYARR, where X = aminosuberic acid) 28 which inhibits Bak, permitting to evaluate the structural and energetic differences that make it and inhibitor rather than an activator of Bak. A mutated form of Bid has also been analysed: BidAib 2 (EDIIRNIARHLAXVGDXMDRSI, where X = Aib), which has increased helicogenic residues, allowing to determine if an increased helical content can improve binding affinity. A non-globular conformation of Bak in complex with Bim 28 and with the inhibitor form of Bim are also studied to determine whether the interaction between Bak and the BH3-only peptide conserves the molecular determinants. These non-globular conformations are called active conformations, as they simulate the structure that are adopted when latch domain is separated from the core domain, in the initial steps during Bak activation. The results reported shed some light into a better understanding of the activation/inhibition mechanism of Bak and could be used for the development of novel activators/inhibitors of pro-apoptotic Bak to be translated into novel therapeutic agents.

Preparation of Bak/X Complexes
Three-dimensional structures of the Bak/Bim complex in its active state and of the Bak/Bid and Bak/Bim(h3Pc) complexes in their inactive state were retrieved from Protein Data Bank (PDB), with PDB codes 5VWV 28 , 2M5B 29 and 5VWZ 28 respectively (Table S1, ESI). When needed, modifications to obtain the wild type (WT) sequence were performed for both Bak, Bim and Bid. The structure of the Bak/Bim complex in Bak inactive conformation was generated using the coordinates of the Bak/Bim(h3Pc) complex (PDB code 5VWZ). 28 The non-natural residue 155 (with residue name 9R1) was modified into its wild-type form: an Isoleucine, and missing atoms from this residue were completed using the LEaP module of the AMBER16 software. 30 Finally, the structure of the Bak/Bim(h3Pc) complex in Bak active conformation was generated from the Bak/Bim complex (PDB code 5VWV) where residue I155 of Bim was changed to the non-natural 9R1 residue.
Parameters for the h3Pc non-natural residue were obtained using distance, angles and dihedrals from glutamate, whilst missing parameters were obtained using the General Amber Force Field (gaff) 31 . Charges were generated with the restrained electrostatic potential (RESP) at HF/6-31G* level, which is the default charge approach applied in Amber protein force fields. 31,32 The Bak/BidAib 2 complex has two special residues, the α-aminoisobutyric acid residues (Aib), which are used due to their strong helicogenic tendency, since they can promote helical folding in synthetic and natural sequences. 33,34 This complex was used to determine whether the use of BH3 mimetic peptides with an increased helical content can improve binding affinity. Moreover, it is important to analyse if these residues have any implication in molecular determinants of the Bak/BH3 interaction. Parametrization of the residue Aib was obtained from Bisetty et al. 35

Molecular dynamics of Bak and Bak/X Complexes
All calculations reported in the present work were performed with the PMEMD (Particle Mesh Ewald Molecular Dynamics) code of the AMBER16 software in its CUDA version 31 using the AMBER ff14SB force field. 36,37 Before performing molecular dynamics (MD) simulations, each system was first relaxed. Systems were soaked in a cubic box, using TIP3P 38 explicit water molecules with periodic boundary conditions. Minimizations were carried out stepwisely: first, positions of the water molecules and ions were optimized using the steepest descent (SD) algorithm 39 up to 5000 cycles of minimization, keeping fixed the rest of the system; in a second step, relaxation of the modelled residues was carried out in two stages, each one of 5000 cycles of SD, keeping fixed backbone positions of the modelled residues and using a decreasing force constant of 5.0 and 0.1 kcal/Å; in the last step, the minimization of the whole system was carried out without any restrictions using 10000 cycles of SD method.
After minimization, systems were heated in steps of 30 K every 20 ps, maintaining all backbone atoms constrained with a force constant of 1.0 kcal/mol·Å. The heating process was performed under the canonical (NVT) ensemble. Once the systems were heated, a 200 ps simulation at constant pressure (NPT ensemble) was performed for density equilibration. The final structures were used as starting points for MD simulations. Trajectories were performed at 300 K by coupling the system to a thermal bath, using the Langevin thermostat. 40 All bonds involving hydrogen atoms were fixed using the SHAKE algorithm, 41 which allowed us to use a time step of 2 fs in all the simulations. Non-bonded interactions were truncated using a cutoff of 9 Å, and long-range interactions were treated with the particle-mesh Ewald summation method. 42 Three MD of 300 ns length were performed on each system to ensure their energy convergence and to enhance conformational sampling, but only the last 50 ns of each trajectory were considered for analysis. Sampling was carried out using multiple molecular dynamics trajectories, which results in a more thorough exploration of the conformational space reducing the chances the system gets trapped in a local minima. 43

Root-Mean Square Deviation (RMSD) and Root-Mean Square Fluctuation (RMSF)
Root-Mean Square Deviation (RMSD) along the simulation time was computed using the cpptraj 44 module from AMBER16 to assess the structural stability of the systems along the Molecular Dynamics simulation. RMSD was computed using the respective experimental structure obtained from LEaP, that is not minimized, as reference structure. Each of the frames of the diverse trajectories were reoriented over the residues from helices α1, α3, α4, α5 and α6. Root-Mean Square Fluctuation (RMSF) for each of the residues of Bak and the diverse peptides were also computed using cpptraj. This provides information of the local conformational flexibility of each residue.

Principal Component Fluctuations
In order to determine and analyse the principal structural variations of the diverse systems studied, we computed the Principal Component Fluctuations of the most relevant secondary structures. For this purpose, we used the Principal Component Analysis (PCA), a multivariate statistical technique systematically applied to reduce the number of dimensions needed to describe protein dynamics through a decomposition process that describe motions from the largest to the smallest spatial scales. In this method, a covariance matrix constructed using the atomic coordinates of the alpha carbons (C of every residue is diagonalized to produce a set of eigenvectors or Principal Components (PC (i) , i=1, N, being N the number of residues of the protein), as well as their corresponding eigenvalues  (i) . Thus, this covariance matrix is a 3N x 3N symmetric matrix. PCA transforms correlated variables into uncorrelated variables and once they are rank ordered, the first principal modes or eigenvectors can be used to characterize large-scale protein motions. In other words, these first modes are enough to define the "essential" space or motions of the protein. 45 The percent contribution of the i th principal component PC (i) , to the structural variance in the dataset is given by , where the summation is performed over all 3N components. The square % = 100 * where the square brackets represent a 3x3 submatrix of the full 3N x 3N matrix and the subscript kk stands for the k th diagonal element of this reduced submatrix. Results of this analysis provide the square summation of the fluctuations of the PCs (PCF) for residues of helices α1, α2, α3 and α5 and each PC fluctuation is computed using the current PC and all the previous PC cumulatively.

Binding Free Energy Calculations
Binding Free energies were computed using the Molecular Mechanics Poisson -Boltzmann Surface Area (MMPBSA) and the Molecular Mechanics Generalized -Born Surface Area (MMGBSA) algorithms, 47 implemented in the AMBER16 package. For both methods, free binding energy can be obtained according to the equation: where ∆H gas is the gas-phase interaction energy calculated by summing the internal energy, the noncovalent van der Waals ( ), and electrostatic ( ) molecular mechanics energies. On the Δ Δ other hand, ∆G solv is computed as the sum of polar ( ) and non-polar terms ( ).

Δ Δ Δ
can be calculated numerically by solving the Poisson-Boltzmann (PB) equation 48 , or the simplified form, the Generalized Born (GB) method 49 , for MMPBSA and MMGBSA algorithms, respectively. In the present work, we used the Onufriev-Bashford-Case (OBC) generalised Born method (igb=2) 50 . Regarding to r , it can be calculated using the following equation: where SASA is the Solvent-Accessible Surface Area, calculated using the LCPO method, 51 and the values for γ and β constants were set to 0,00542 kcal/mol·Å 2 and 0,92 kcal/mol for MMPBSA 52 and 0,0072 kcal/mol·Å 2 and 0 kcal/mol for MMGBSA 53 . The entropic term was calculated within the harmonic entropy approximation using the normal-mode analysis with default values as implemented in the MMPBSA.py program 54 . Because of the high computational requirements, the entropic contribution was averaged over 80 snapshots extracted at a time interval of 250 ps from the last 20 ns of each trajectory.
In the present work, binding free energies of the Bak/Bid, Bak/Bim, Bak/BidAib 2 and Bak/Bim(h3Pc) complexes were computed from three MD trajectories using MMPBSA and MMGBSA and the energy values were averaged over 50 ns after ensuring that each of these trajectories was stable.

Analysis of Hydrogen Bonds
Hydrogen bonds between the Bak protein and the BH3 peptide ligands were also analysed using the cpptraj module of AMBER16. The presence of these bonds was calculated over the last 50 ns of simulation time, and only hydrogen bonds with occupancies greater or equal to 60% were considered. Criteria used for hydrogen bonding was that the distance between the hydrogen-bond acceptor and the hydrogen-bond donor was lower than 3.0 Å and the angle between Acceptor -H -Donor more than 135 .°C lusterization Similar structures for every complex were grouped into 10 different clusters using the cpptraj module of AMBER16. This process was performed using the average link algorithm 56 and the backbone Cα root-mean-square deviation of the superposed residues, the α3 and α4 helices involved in the interaction between Bak and the BH3 peptides, was used as a measure for the distance between two conformations

Results
As mentioned above, the aim of this study was to get a better understanding of the characteristics that allow the interaction between Bak and specific BH3 peptides (Bid, Bim, BidAib 2 and Bim(h3Pc)) and to analyse the structural and energy alterations that are taking place in all the systems. We have taken special interest in analysing the differences between Bak/Bim and Bak/Bim(h3Pc) complexes, being the former an activator of Bak and the latter an inhibitor, aiming at determining structural and energy differences that produce the change in activity. 28 We also wanted to analyse the differences between Bak/Bim and Bak/Bim(h3Pc) complexes on their active conformations, and determine if those differences that occur in the globular structures can be extrapolated in the active conformations, or in contrast, a different behaviour is observed. When analysing these results, structural and energy stability from different complexes were considered. Structural stability of each complex was calculated using the RMSD method and the analysis of the fluctuation has been performed with RMSF methods and with the analysis of the fluctuation of the Principal Components from the system. Binding free energy from these complexes was calculated using the MMPBSA and MMGBSA approaches. Moreover, detailed binding free energies from Bak residues were calculated using the free energy decomposition method.

Root-Mean Square Deviation and Root-Mean Square Fluctuation
RMSD was calculated along the trajectory to establish if the different MD simulations performed for every complex converge and are stable in structural terms. Residues superposed to perform these analyses were residues from helices α3, α4 and α5, involved in the binding of BH3-only peptides, and those from α1 and α6, which are relevant in the conformational changes that Bak experiences when it is activated. In the active conformations, α6 was not considered since it experiences an increase in its mobility due to its exposed conformation.
As it can be seen in Figure 2, the RMSD of the diverse Bak/BH3 peptide complexes experience an increase at the beginning of the trajectory, but after 50ns, all systems are structurally stabilized. When referring to the active conformations, they show a higher RMSD, as the latch domain is not  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 interacting with the core domain, and the globular structure moves to fill the space left by this latch domain. RMSD for the Bak/Bid complex is slightly higher than that for the other globular conformations. However, these differences are not as significant as those observed between active and inactive conformations.
An analysis of the root-mean-square fluctuation (RMSF) per residue for the diverse complexes (Figure 3), allows to conclude that all four globular complexes exhibit similar RMSF distributions, showing a rigid behaviour in almost all the structure, with the α1 -α2 loop showing the largest fluctuations. On the other hand, the active conformations show much higher fluctuation, due to the movement of the core domain, which is trying to fill the space left from the latch domain.

Principal Component Fluctuation
The analysis of RMSF failed to identify differences between activators and the inhibitor in structural terms. Thus, we decided to perform PCA, a powerful method to reveal the principal variations in structure, studying the fluctuation of different residues (PCF). Figure 4 shows the computed square summation of the atomic fluctuations for helices α1, α2, α3 and α5. As can be seen, no differences between Bak/Bim and Bak/Bim(h3Pc) can be observed for residues of α2 and α5 ( Figures 4B and 4D respectively). However, for the last residues of helix α1 (residues from 42 to 46) ( Figure 4A), the activator exhibits larger fluctuations in comparison to the inhibitor. This difference is due to the interaction between residues 9R1 155 from Bim(h3Pc) and R42 from Bak, which is not found in Bim wild type (WT). This difference in fluctuations may modify the structure of helix α1 as well as the loop between α1 and α2, which is in close contact to α6, that is involved in the activation process of Bak. As it can be seen, addition of more Principal Components does not improve the explanation of the fluctuations of the system. Moreover, it was recently identified an activation site in this loop triggered by antibodies. 24 Thus, if the inhibitor is decreasing the fluctuation of this region, it could be affecting conformational changes that are necessary to produce a full activation of Bak.
Other differences between the activator and the inhibitor can be seen in diverse residues of helix α3 ( Figure 4C). In the complex with the activator, residues 87-92 experience large fluctuations, whereas in the complex with the inhibitor these are found between residues 90-95. This difference is due to the interaction between 9R1 155 from the inhibitor and N86, reducing the fluctuation at the beginning of the helix. When Bak is on its inactive form, helix α3 is stabilized by the interaction between D90 and R42. In the complex of Bak and the inhibitor, this interaction is being replaced by the interaction of 9R1 155 of Bim(h3Pc) with R42, destabilizing helix α3.
The analysis of the fluctuation of the Principal Components of Bak/Bim and Bak/Bim(h3Pc) in the active conformations show significant differences in regard to the inactive conformations ( Figure  5). Large fluctuations can be seen in the whole structure as was observed in the RMSF analysis. Fluctuation of residues of helix α2 shows no differences between Bak/Bim and Bak/Bim(h3Pc), as it can be seen in Figure 5B, the same as observed in the inactive conformations. However, a large fluctuation of residues in helix α5 is found in the Bak/Bim complex compared to the Bak/Bim(h3Pc) one ( Figure 5D). This increase of residue fluctuations is a consequence of the separation of the latch domain from the core domain, since the absence of α6 produces an increase in the fluctuation to occupy the free space. However, in this active conformation, the inhibitor form of Bim is fixing α5 helix through its interaction with R137. This decreased fluctuation may be involved in the suppression of later conformational changes needed for complete Bak activation. Similar behaviour is experienced by residues from helix α1 of the Bak/Bim complex (see Figure  5A) in comparison to the Bak/Bim(h3Pc) complex, where larger fluctuations are observed in all the residues and not only in the last as observed in the inactive conformations. In this case, in addition to the interaction formed between Bak R42 and 9R1 155 from Bim(h3Pc), the separation of the latch domain contributes to the increase in the fluctuation of helix α1 residues. As happened for helix α5, this decreased fluctuation in the Bak complex with Bim(h3Pc) in comparison to the complex with the WT form of Bim could affect conformational changes that take place later to fully activate Bak.
As shown in Figure 5C, residues of helix α3 fluctuate in both Bak/Bim and Bak/Bim(h3Pc) complexes. However, the analysis of the average appearance frequency of hydrogen bonds between Bak residues and 9R1 residue of Bim(h3Pc) in the inactive and active complexes ( Table 1), reveals that, in the active conformation, interaction between residue 9R1 155 from Bim(h3Pc) and N86 does not occur with the same frequency as in the case of inactive conformation. Thus, the first residues from helix α3 show a higher fluctuation in the inactive conformation than in the active from Bak/Bim(h3Pc) complex. Focusing on the last residues of helix α3, fluctuation is slightly higher in Bak/Bim complex than in Bak/Bim(h3Pc). This explanation could be contradictory with what it is seen in the inactive conformation. However, in the active conformation from Bak/Bim complex, helix α1 is not stabilizing helix α3, producing the increased fluctuation observed.
This difference is due to the interaction between 9R1 155 from the inhibitor and N86, reducing the fluctuation at the beginning of the helix. When Bak is on its inactive form, helix α3 is stabilized by the interaction between D90 and R42. In the complex of Bak and the inhibitor, this interaction is being replaced by the interaction of 9R1 155 of Bim(h3Pc) with R42, destabilizing helix α3.

Binding Free Energy of the Complexes
Binding free energy was computed using MMPBSA and MMGBSA approaches for all the complexes along the trajectory to determine the convergence of the different trajectories, their stability in terms of binding energy and to carry out a more exhaustive analysis using the last nanoseconds of the simulations. As it can be seen in Figures S1 and S2, the diverse complexes have different ∆G binding along the trajectory. Complexes between Bak and BH3 peptides have a similar ∆G binding , with values between -90 and -100 kcal/mol. However, the Bak/Bim(h3Pc) complex exhibits a more negative ∆G binding , around -140 kcal/mol. This increased ∆G binding suggests a higher stability of the complex, preventing the dissociation of the BH3 peptide from Bak. Table 2 lists a detailed decomposition of the binding free energy computed using the MMPBSA and the MMGBSA approaches for the diverse complexes using the last 50 ns of the corresponding MD trajectories. As it can be seen, van der Waals contribution (∆H vdW ) is notably important to favour binding. It is important to mention that total electrostatic contribution (∆H ee + ∆G polar ), regardless of the approach used, is unfavourable for the binding of all complexes, which underlines the importance of the ∆H vdW contribution. On the other hand, the entropic contribution to the binding energy is similar for all the complexes studied. However, its contribution in the case of the inhibitor system, Bak/Bim(h3Pc), is higher than the corresponding to the Bak/Bim activators, for both the active and the inactive complexes. This suggests that the introduction of non-natural amino acids is not important in terms of the entropic contribution to the binding energy.
The most significant result of these calculations is the higher binding free energy of the Bak/Bim(h3Pc) complex in comparison to the rest. This difference is due to the presence of the 9R1 residue in Bim(h3Pc), which contributes to increase substantially the electrostatic term, which results in this increased ∆G Total . This phenomenon can be also observed in the active conformations, in which Bim(h3Pc) total free binding energy is 20 kcal/mol higher than that of Bim. In the case of Bak/Bid and Bak/Bim, the average binding free energy computed during the last 50 ns are similar. Similar values for ∆G TOTAL are obtained in the case of the Bak/BidAib 2 system, despite their increased helicogenity due to the presence of the Aib residues.
When comparing active and inactive conformations without the entropic contribution, no significant differences can be observed in the Bim complexes. In the case of the inhibitor Bim(h3Pc) systems, the inactive conformation appears more stable. However, these trends are different when the entropic contribution is included. In this case, when the inhibitor is bound, the two conformations are similar; in contrast, for the activator there is a small difference between the two conformations. Moreover, as shown in Figures S3 and S4, not only the Bak/Bim complex exhibits a similar energy profile during all molecular dynamics simulations in both conformations, but also Bak/Bim(h3Pc) exhibits a similar profile. Thus, despite the large difference in structural terms observed in the analysis of RMSD, in energy terms these differences are more subtle, suggesting that the separation of the latch domain from the core domain implies a small difference in the interaction between Bak and an activator or an inhibitor and in this process the entropic changes can play a subtle role.
It is generally accepted that in the complexes formed between antiapoptotic and proapoptotic proteins with BH3 peptides, the most relevant residues in terms of energy involved in the interaction are hydrophobic. This observation is also seen in Figure 6 and Table S1, where for Bak bound to direct activator BH3 peptides, the residues with higher energy contribution are mainly hydrophobic. The exception is R127, whose contribution is mainly electrostatic due to it is involved in a hydrogen bond with BH3 peptides (presence of this hydrogen bond, higher than 70% in all systems analysed).This hydrogen bond not only occurs in the interaction between pro-apoptotic and BH3-only proteins, but also in the interaction between anti-apoptotic and BH3-only proteins. 57 Comparison of the diverse free energy contributions per residue of the diverse complexes with activators with those of the complex Bak/Bim(h3Pc), shows that some residues exhibit an important contribution as a consequence of the I155X (X = suberic acid) mutation of Bim. This mutation produces an approximation to residues of helices α1, α3 and α5 that were not reached by the Isoleucine. These residues are R42, N86 and R137, which are establishing a hydrogen bond with the 9R1 155 residue from Bim(h3Pc), as it can be seen in Figure 7. The frequency of the hydrogen bonds involving this residue can be seen in Table 1.
This per residue free energy contribution is similar when comparing inactive and active conformations, indicating that interactions that take place in the globular structure also take place in the active conformations. Thus, the study of the interaction between BH3-only peptides with proapoptotic proteins can be performed in the globular form, as it has been performed until nowadays, without modifying the interaction pattern.

Pharmacophore Definition
In order to identify the stereochemical features that characterize the binding of different BH3 peptides, the most populated representative for each complex obtained from the clusterization of all three MD simulations was used. Residues from the BH3-only peptides with the higher energy contribution to the interaction have been used to design the pharmacophore. As it can be observed  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 in Figure 8, the pharmacophore is very similar in all complexes. BH3 peptides bind to the canonical hydrophobic groove (α3, α4 and α5) of Bak primarily through hydrophobic features. A conserved Hydrogen Bond interaction is present between Bak-R127 and BH3-D95/D157 (Bid and Bim respectively). The mutated form of Bim, Bim(h3Pc), differs in the third hydrophobic feature, where appears a Hydrogen Bond acceptor feature in the canonical hydrophobic groove, in a cavity from this groove.
Analysis of the interaction between Bak and BH3-only proteins has always been performed in the globular and inactive form of Bak. Thus, we wanted to analyse if there is any difference in the most relevant features in the interaction when comparing this inactive form of Bak to its active conformation. Results of distances and angles between these molecular features are shown in Figures  S5 and S6 and Tables S2 and S3. In these results it is clearly seen that there are no relevant differences in the most important residues in the interaction between Bak/Bim and Bak/Bim(h3Pc) either in the active or the inactive conformations. Thus, it confirms what was expected: if we want to find new molecules that bind Bak, the study can be developed using the active or the inactive conformation indistinctly.

Conclusions
Understanding the mechanism by which pro-apoptotic Bak is activated is a promising strategy for designing an effective cancer therapy, as it develops a key role in the crucial process of apoptosis. It is widely known that Bak remains inactivated in some diseases. Thus, the study of Bak structure and interactions present in Bak protein complexes with BH3 peptides should give some perspective for designing small molecules that can bind to Bak as other BH3-only proteins do. Moreover, inhibition of Bak to prevent mitochondrial disruption may also be an effective strategy to inhibit cell death.
The analysis of the structural differences between Bak/Bim and Bak/Bim(h3Pc) showed that, for globular or inactive structures, differences are present in helices α1 and α3, due to the interaction formed between 9R1 residue of Bim(h3Pc) with R42 and N86 respectively. Fluctuation shift in the last residues of α1 could explain the fact that conformational changes that take place during Bak activation are being avoided. This conformational change is the dissociation of α1 helix from α6, which is needed not only to separate the latch domain from the core domain, but also to form the apoptotic pore. On the other hand, the different fluctuation between Bak/Bim and Bak/Bim(h3Pc) in α3 helix can be explained by the interaction between 9R1 with N86, but also with the fact that the interaction between 9R1 and R42 is competing with an intramolecular hydrogen bond between this arginine and D90, from α3 helix, increasing the fluctuation of this helix in the inhibitor as seen in Figure 4C.
As the activation mechanism of Bak, and more accurately, the conformational changes that Bak needs to experience to insert in the MOM are not fully understood, the analysis of an active conformation of Bak has been performed. These active conformations resemble a posterior step of Bak activation, where latch domain, helices α6 to α8, has been separated from the core domain, helices α2 to α5. As all previous analyses of Bak had been performed using the globular structure, we wanted to analyse if there are any differences in the interaction between active and inactive conformations. In terms of structural differences, fluctuation is increased in overall structure, as the separation of latch domain leaves an empty space, filled by the helices from the core domain. However, this increase in the RMSF is not observed in helices α1 and α5, in which the interactions with 9R1 residue stabilize and decrease the movement of these helices. This reduction in the fluctuation is showing that, despite the latch domain is already separated from the core domain, posterior conformational changes in further Bak activation may be inhibited.
Moreover, we wanted to analyse Bak/Bim(h3Pc) complex on its active conformation to determine the inactivation mechanism of this inhibitor. It has become clear that the inhibitor is causing a difficulty in performing Bak conformational changes during activation. However, there is a possibility that this inhibitor is also preventing Bak dimerization. Energetic analysis of the interaction between Bak and BH3-only peptides has been carried out and it has been observed that the complex with the inhibitor has a better free binding energy than with any other activators analysed. We have hypothesized that this increase in the interaction indicates a higher difficulty in dissociate Bak from BH3-only peptide, preventing the association of Bak-Bak by their BH3 domains. Nevertheless, to assure this hypothesis, an analysis of the Bak-Bak interaction must be performed.
Predictions of binding free energy values are a key issue for determining the binding pattern of a given molecule. These predictions are a complex issue, especially when structural data is obtained from molecular dynamics simulations. Hence, two different approaches were used, MMPBSA and MMGBSA, to have a proper performance when estimating binding free energy. Despite the expected fluctuations observed in these predictions, both approaches and time intervals (last 50 ns of each MD simulation) analysed lead to stability and equivalent results (Figures S1-S4).
Binding of BH3-only peptides to Bak protein is mainly conducted by hydrophobic interactions. These hydrophobic features exist in BH3-only peptides conserved positions, confirming that BH3 domain is constituted of conserved hydrophobic residues in conserved positions. Residues from these positions are: for Bid and BidAib 2 , I86, L90, V93 and M97; and for Bim and Bim(h3Pc), I148, L152, I155 (9R1 for Bim(h3Pc)) and F159. Due to this major presence of hydrophobic features in the interaction between Bak and BH3-only peptides, finding small selective molecules to activate Bak may signify an increased difficulty.
It is true that hydrophobic interactions govern substantially the interaction between Bak and BH3-only peptides. However, a hydrogen bond is present in all molecular dynamic simulations from all complexes. This hydrogen bond, with R127 of Bak, is key in order to find more selective molecules to activate Bak. Thus, it can be affirmed that although the binding of BH3-only proteins to Bak is mainly hydrophobic, electrostatic interactions might also contribute to complex formation and increasing the affinity of these proteins. Taking that into account, D157 in the case of Bim and D95 for Bid, has been found critical for the interaction with R127 of Bak in all complexes. However, this interaction is not only present in pro-apoptotic proteins, but it is also conserved in the interaction between BH3-only proteins and antiapoptotic proteins. Therefore, this interaction is not specific of pro-apoptotic proteins.
Analysis of the energy contribution by each residue was performed to determine the differences in the interaction between Bak/Bim and Bak/Bim(h3Pc). While in Bak/Bim complex interactions from the canonical hydrophobic groove are mainly hydrophobic, except for the hydrogen bond between R127 and D157 of Bak and Bim respectively, in Bak/Bim(h3Pc) complex, this hydrogen bond exists as well as the hydrophobic interactions. However, what increases affinity between Bak and Bim and produces the already commented change of activity, inactivating Bak, is a hydrogen bond formed between R42 or R137 from Bak and 9R1 residue from Bim(h3Pc).
Interaction between 9R1 and R42, R137 or N86 is really important for the activity switch of Bak. Structural evidences in terms of fluctuation have demonstrated this importance. These evidences are supported by the energetic results, which show an increased affinity for the inhibitor form. Thus, this interaction helps decreasing the movement of these residues, which results in a more "fixed" conformation.
All these changes that Bak suffers when is interacting with Bim(h3Pc) produce the previous mentioned inactivation of Bak. As Bak structure and energetic terms are modified, two possibilities arise that could inactivate Bak. On the one hand, as the movement of Bak is being altered, it could be that what the inhibitor is doing is preventing conformational changes that take place during Bak activation, which could be either the dissociation of the α1 helix 19 or the separation of core and latch domain. As differences between Bak/Bim and Bak/Bim(h3Pc) also occur in the active conformations, it could be that posterior conformational changes to fully activate Bak are also altered by the inhibitor. On the other hand, as the affinity between Bak/Bim(h3Pc) is considerably higher than the one from Bak/Bim complex, it is possible to assert that the dissociation of the inhibitor will be more difficult, and the transient binding of BH3-only peptides would be more permanent, obstructing the dimerization of Bak. In both possibilities the ability to permeabilize the MOM and to release cytochrome c from the mitochondria is blocked, so apoptosis does not occur.
Pharmacophore features in the interaction between Bak and BH3-only peptides are nearly the same when comparing inactive and active conformations. Only some differences in distances and angles that involve residue 155, which is the one that varies, can be seen. However, these differences are not significant, as it can be seen in Figures S5 and S6. Thus, it is possible to conclude that the analysis of the interaction between Bak and BH3-only peptides can be performed in either the inactive conformations, the way it has been performed so far, or the active conformation, with the latch domain separated from the core domain, a more similar conformation when Bak is inserted in the MOM.
Development of agents that are directly able to activate or inhibit Bak are important likewise: further investigation in the activation of Bak by BH3 peptides is crucial to improve cancer therapy; and inhibiting Bak might open many prospects to a better management of degenerative disorders.
Increasing knowledge about the activation mechanisms of pro-apoptotic proteins, and more accurately, of Bak, is key to understand how cells develop the crucial process of apoptosis. Further analysis and studies on Bak structure and on its interaction with BH3 peptides are needed, as well as the study of Bak dimerization and oligomerization. Understanding how Bak is modified during its activation and how it dimerizes and oligomerizes is essential in order to find new molecules that simulate all these modifications and activate Bak in those cases that it is inactivated.      Figure 1: Crystal structure of Bak in complex with a BH3-only peptide. Core domain is formed by helices α2, α3, α4 and α5, and latch domain is formed by α6, α7 and α8. Full length pro-apoptotic Bak has 9 α-helices but α9, which is the transmembrane helix, has not been crystalized. (A) Bak on its soluble or globular form in complex with a BH3-only peptide. (B) Bak in complex with a BH3-only peptide in a conformation of the early activation steps, with the two domains, core and latch, separated.      1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58 59 60