BIM Binding Remotely Regulates BAX Activation - ACS Publications

Dec 26, 2017 - mostly by minimizing the free energy difference between the active (α9-released) and inactive (α9-bound) states. ..... 80°) energy i...
0 downloads 9 Views 2MB Size
Subscriber access provided by the University of Exeter

Article

BIM Binding Remotely Regulates BAX Activation: Insights from the Free Energy Landscapes Souvik Sinha, Atanu Maity, and Shubhra Ghosh Dastidar J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00628 • Publication Date (Web): 26 Dec 2017 Downloaded from http://pubs.acs.org on January 2, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40 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

Journal of Chemical Information and Modeling

BIM Binding Remotely Regulates BAX Activation: Insights from the Free Energy Landscapes Souvik Sinha1, Atanu Maity1, Shubhra Ghosh Dastidar1,* 1

Bioinformatics Centre, Bose Institute, P-1/12 CIT Scheme VII M, Kolkata 700054

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

Abstract: Activation of the pro-apoptotic BAX protein, a BCl-2 family member, is known to trigger the apoptosis by forming pores in the mitochondrial outer membrane (MOM). While in the cytosol, release of its transmembrane C-terminal helix (called α9 helix) from a well-characterized binding pocket (BC groove) and subsequent permeabilization of MOM is understood to be the initiating events of the activation. In concern of what initiates the BAX activation, so far one plausible suggestion has been that the transient attachment of BH3-only peptide on a distal site from BC groove triggers the activation process. Yet how this pivotal step displaces α9 from the BC groove has remained poorly understood. Using a combination of standard molecular dynamics and enhanced sampling methods, the energy landscape of BIM (BH3-only peptide) induced BAX activation has been computed and molecular origin of those events is hereby reported in atomistic detail. The simulated transition pathway of α9 release reveals that BIM subdues the energetic cost of the process by reducing the activation energy barrier to some extent but mostly by minimizing the free energy difference between the active (α9-released) and inactive (α9bound) states. Interestingly, the flexibility of the α9 helix itself plays a decisive role in this mechanism. The impact of BIM encounter at the distal site is found to propagate to the α9 (BC groove bound) mostly through conserved pathways of residue level interactions. Overall, the thermodynamic basis of the ‘hit-and-run’ mechanism for activation of BCL-2 family is presented reconciling the available biochemical observations.

2 ACS Paragon Plus Environment

Page 2 of 40

Page 3 of 40 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

Journal of Chemical Information and Modeling

Introduction: Mitochondria, the energy factory of a cell, harbor the electron transporter hemeprotein Cytochrome C (Cyt C) whose release from mitochondria triggers the apoptotic death of the cells. The BCL-2 family regulates the mitochondrial pathway of apoptosis and this is controlled by the interplay between pro-apoptotic (e.g. BCL-XL, BCL-W, BCL-2 etc), anti-apoptotic (e.g. BAX, BAX) and BH3-only (e.g. BID, BIM etc.) proteins of the family. The pro-apoptotic sub-group of BCl-2 family facilitate the process of releasing Cyt C from mitochondria by forming pores in the mitochondrial outer membrane (MOM) when required; these proteins otherwise get antagonized by the anti-apoptotic counterpart of the same family. Although these subgroups have complementary functions, a substantial portion of their structures share 3 to 4 homology domains (indexed as BH1, 2, 3, 4 domains), whereas activator proteins share only the BH3 homology domain.1-3 The BH3-only proteins bring the death signal from different cell stresses to the apoptotic machinery, either get neutralized by anti-apoptotic members or activate the proapoptotic proteins, depending on the cellular requirement.4, 5 Substantial evidence is available for the activation of pro-apoptotic BAX/BAK by direct binding of BH3-only proteins,5-11 which further proceeds to trigger mitochondrial outer membrane permeabilization (MOMP). Though inactive BAK resides at MOM, inactive BAX exists as a monomer in the cytosol and upon activation by BH3-activator proteins, it integrates into MOM to induce MOMP.5, Monomeric BAX consists of nine α-helices [Figure 1(A)]

17

12-16

(numbered as α1, α2, α3 etc. up to

α9 from N to C terminal). The range of residues correspond to each helix is presented in TableS1. The C-terminal helix (α9), which appears to be a transmembrane domain, also fits into the well-characterized hydrophobic BC-groove (BCG) (a pocket known to accommodate either BH3 peptides as ligands or the α9 as pseudo-substrate) in the inactive form of BAX, but gets exposed 3 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

Page 4 of 40

to MOM for anchoring upon activation7, 18, 19 [Figure 1(B)]. The release of α9 from the BCG is the primary step of BAX activation which then proceeds to MOMP followed by further conformational alterations in the cytosolic globular domain. Such alteration enables activated the activated monomers to oligomerize to form pores in the MOM through which Cyt C and several other apoptogenic factors (Smac/Diablo, Omi/HtrA2 etc.) get released.15, 16, 18, 20, 21 What activates the BAX-like protein and how it does so are the paramount questions of apoptosis research; the answers are believed to be ‘the Holy Grail of apoptosis process’.3 It is well documented that the initial steps of activation of BCL-2 family members involve the release of α9 from the BCG to make room for an incoming BH3-only protein which then binds there to trigger the activation.22, 23 Even the BAK activation is reported to be facilitated by the binding of certain BH3-only peptide to the canonical surface groove i.e. the BC groove.24 In contrast to this well-accepted understanding, several experiments with full length BAX protein7,

19, 25

have

hypothesized that the release of α9 in cytosolic BAX monomer is eased by an allosteric push originated from the rear side pocket (defined by α1 and α6, see Figure 1B) of the canonical BC groove upon ‘transient’ encounter with BH3-only proteins (e.g. the BIM peptides); the exact molecular mechanism of such distal activation is yet to be clearly understood. However, with α9truncated BAX, Czabotar et al.26 provided evidence of activator BH3-only protein binding at the hydrophobic BC-groove (in ‘membrane-embedded’ population) and causing dislodgement of the ‘latch’ (α6 - α8) domain from the ‘core’ (α1- α5) domain for further dimerization. Replacement of the activator BH3-only proteins by endogenous BH3-domain of another activated monomer is important to form ’BH3-in-groove’ dimers27-31 of activated pro-apoptotic BAX or BAK proteins. This dimeric complex further accumulates more protein units to construct the pore-forming oligomer in MOM. 4 ACS Paragon Plus Environment

Page 5 of 40 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

Journal of Chemical Information and Modeling

The state-of-art experiments are yet to reveal the molecular origin of BAX activation in atomistic details. There, computational modeling from reported experimental data (e.g. solution structures, biochemical reports etc.) together with molecular dynamics (MD) offers a quality opportunity for such detail investigation. The MD simulation can further be used to extract the multidimensional energy landscapes associated with complex biological events.32-36 As experiments suggested multiple layers of hypothesis on the BAX activation process, it would be easier to catalog them with the knowledge of the energy surface correspond to each hypothesis and rationalizing each step of activation with shreds of evidence. Therefore using computational approaches, the present investigation has estimated the free energy landscapes of BAX activation in terms of α9 release from the BCG (considering this as the first major structural modification towards the activation). Furthermore, energy landscapes in presence and absence of BIM peptide (BH3-only) at the rear pocket of a cytosolic BAX monomer are compared and the associated structural consequences are characterized in atomistic details. This study also rationalized the thermodynamic requirement of ‘transient’ binding of the BH3-only peptide for the sake of BAX activation.

Methods: Structural Modeling: The solution structure of cytosolic BAX monomer (PDB ID: 1F16)17 was used as the starting model of apo BAX (ABX) [Figure 1(A)]. For the holo BAX:BIM complex (BBX) [Figure 1(B)], BIM BH3-peptide coordinates were taken from 2K7W25 and BAX monomer coordinates were modeled using 1F16 except its loop joining α1 and α2 (residue 3854). The loop joining α1 and α2 was modeled using the 2K7W as a template to maintain the open conformation of the loop in the complex (free outside the α1/α6 pocket) as suggested by Gavathiotis et al.25. Moreover, the α1-α2 loop starts at Met38 in 2K7W (BAX:BIM complex), 5 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

and at Gly36 in 1F16 (apo BAX monomer). This difference was kept in mind during modeling of the loop in respective starting structures. Apart from this ~20 residue α1-α2 loop, there is a 15 residue long N-terminal flanking end (residue 1-15) which caused larger RMSD (~6-8 Å) in previous full-length BAX simulations for tens of nanoseconds37, 38 to even 1 millisecond39, yet part of the N-terminal flanking end (Gly12-Ala24) is important for BAX cellular distribution.26, 40

Considering both the difficulty of convergence and importance of the N-terminal for cellular

distribution, only first 9 residues from the N-terminal end of BAX monomer were truncated for modeling of ABX and BBX systems to get sufficient details from reasonable simulation length (hundreds of nanoseconds). All the starting structures were prepared using PyMOL software.41 The models were then refined using molecular dynamics (MD) simulations. MD Simulation: In each system the macromolecules were solvated by a cubic box of TIP3P42 water models, ensuring at least 10Å thickness of water layer everywhere. After neutralizing the overall charge of the system by adding ions, excess K+ and Cl- ions were added to top-up the salt concentration to 0.15(M). Modeling and energy minimizations were performed using CHARMM biomolecular simulation program.43 The molecular dynamics of each system was run for 200 ns using NAMD 2.1144 and CHARMM all3645 parameter set (Table-S2). For all these simulations, short-range non-bonded interactions were calculated using cut-off at 14 Å whereas the longrange electrostatic interactions were computed using Particle Mesh Ewald method.46 Temperature was kept constant at 300 K using Langevin dynamics, with a damping coefficient of 1 ps-1. The pressure was also kept constant at 1atm using Langevin piston method47 with a coupling constant of τp = 0.5 ps. The integration time step was set to 2 fs after constraining all bonds involving hydrogen using SHAKE algorithm.48 It is reasonable to assume that the 200 ns MD simulation was sufficient to optimize the models (RMSDs are discussed later) and therefore 6 ACS Paragon Plus Environment

Page 6 of 40

Page 7 of 40 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

Journal of Chemical Information and Modeling

the last frame of each trajectory was used as the starting point of the following sets of Adaptive Biasing Force (ABF) based MD simulations (ABF-MD) to compute the Potential of Mean Force (PMF) profiles of the events of α9 release from the BCG. ABF-MD and PMF calculation: ABF-MD does sampling within predefined bins along the reaction coordinates of interest. The Adaptive Biasing Force (ABF) algorithm is implemented in the Collective Variable module49 of the NAMD. The principle of ABF is to remove the free energy barriers adoptively along a reaction coordinate to accelerate transitions between states without affecting the very nature of an equilibration dynamics, even the random fluctuation force.50 During the production dynamics, an external biasing force equivalent to the running time average of the instantaneous force is applied to cancel out the current estimate of the average force. As the simulation converges over time, the biased average force converges around zero and system experiences free diffusion along defined Reaction Coordinates (RC – these are defined as the most relevant degrees of freedom that can best describe any required conformational change). The following equation presents the working principle   =      ξ(X)-z  

(1)

where  is the free energy along z-coordinate,  is potential which changes as   −  t[ ] and  t is an update that eventually converge with free energy . The Dirac delta function ξ(X)-z   is defined by the subset  ,   =  where  is the given reaction coordinate. Here two RCs were defined [Figure 1(C,D)] for the PMF calculation. (i) Angle of opening (Ѳ) (ANG) between α9 and the hydrophobic core (Cα of Ile66, Thr172 and Gly179 were used for this RC). As the angle increases, α9 leaves the BC groove. This RC was varied from 80˚ to 140˚ 7 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

where 80˚-100˚ range was recognized as the BC-groove residence. This RC is stratified into 6 windows of length 10˚ each as the time of convergence is proportional to the square of the length of the transition path.51 (ii) End-to-end distance (E2ED) of α9 (distance between backbone C of Trp170 and Trp188) was described as another RC. This second RC was chosen for two reasons: [a] defining more than one RC reduces the possibility of trapping under non-ergodic scenario50 and [b] effect of BIM binding on the flexibility of α9 can be probed by this choice as reported in previous studies.17, 39 This RC was varied in the range of 25 Å - 31 Å and stratified into 6 bins of length 1 Å. For each of the two systems (ABX and BBX), ABF windows were simulated for 40ns using Multiple Walker strategy

52

and in each run 500 MD steps were used prior to

applying biases. Each window was divided into 4 replicas where replicas were considered as different “walkers” those explore the phase space independently and shared the information regarding sampled instantaneous forces of the explored regime with each other periodically. This is a useful strategy to accelerate sampling along orthogonal slower degrees of freedom and get rid of metastable trapping.50 Each replica has been simulated independently for 10ns and sampled forces were shared among other walkers with a fixed interval of 4ps. So total 480ns of ABF run was used to construct PMFs. Two additional restraints have been carefully introduced to reduce unwanted interference along other degrees of freedom that can contaminate the desired the energy landscape. A dihedral restraint of 5 kcal/Å2mol was used at the groove: α9 interface to reduce horizontal sweeping of the α9 when it is out of the groove and a harmonic restraint of 20 kcal/mol was used at the BCG backbone Cα to reduce the extensive flexibility of the α3-α4 loop region. All the starting structures along angle RC were generated using CHARMM program. Later 2D PMF is integrated along angle RC to get 1D nature of the landscape using the following equation53:

8 ACS Paragon Plus Environment

Page 8 of 40

Page 9 of 40 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

Journal of Chemical Information and Modeling

e = e

$%&'(, )*

 , ! "  #

dy

$%&'( , )* + dy #

(2)

Where β = (1/kBT) ; kB is Boltzmann Constant, T is absolute temperature, W(x,y) is the 2D PMF for the range of (x,y), W(xc,yc) is value of the PMF at any arbitrary point (xc,yc), W(xc,y) is value of the PMF for an arbitrary x-coordinate but throughout the range of y-axis. Principal component analysis: The principal component analysis (PCA) was performed on Cα atoms on the 200ns trajectories of ABX and BBX using ProDy interface.54 PCA uses that actual movement of the atoms recorded in the trajectory files of the MD simulation and determines the linearly independent components (called principal component or PC) of internal dynamics of the macromolecules

55

. It indicates their relative contribution (% contribution) to the overall

dynamics (total variance) of the system and thus it significantly simplifies the complex simulation trajectory data for understanding.56-59 For an analysis of N number of atoms, only 3N6 modes construct the internal motion (leaving 3 translational and 3 rotational modes) of the molecule and it is generally found that only 3-7 modes determine ~90% of the total variance of the system; therefore those modes constitute the ‘essential dynamics’ of the system. These PCs are ranked [PC1, PC2, to PC(3N-6)] according to their decreasing order of contribution to the total variance of the system. Community analysis and shortest Path calculation: A network is a set of nodes connected by edges where an edge is considered between a pair of nodes (Cα of residues) only when 2 heavy atoms from the residues are within 4.5Å of each other for at least 75% of the length of the trajectory. Later, Girvan-Newman algorithm

60

was used to partition the contact networks into

correlated communities where the nodes (amino acids) within the same community interacts stronger than the nodes of other communities. Shortest paths among specific nodes are also 9 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

calculated by employing Floyd-Warshall algorithm. All the calculation of Network analysis was performed using the Network View plugin 61 in VMD. 62

Results: The trajectories of the standard MD simulations, i.e. MDABX and MDBBX, revealed that the modeled structures were stable as they did not go through any major conformational change and also their respective RMSDs were found to be within the 1.4Å – 2.5Å range (Supporting information Figure S1). The comparison of RMSFs of these two systems did not show any remarkable difference; the only difference lies around C-terminal of α9 and α5 - α6 junction, yet it is not of much significance as these regions are unstructured (Figure S2). In both the systems, α9 was found to be packed within the BCG throughout the 200 ns simulation run and therefore the last structures from the 200 ns trajectories would be the reasonable starting point for the computation of activation PMF profiles. All the observations on the free energy profiles (ABFMD simulations) and equilibrium MD trajectories are discussed in the following sections.

Free Energy Profile: The free energy profiles of α9 release from the BC groove (BCG) as a function of α9 helical length (E2ED) [Figure 1(D)] and the angle between α9 and hydrophobic core (ANG) [Figure 1(C)] are presented in Figure 2. Details of the test of convergence63-65 of PMF and validation of the corresponding transition states from committor analysis are provided in the Supporting Information (SI). As shown in Figure 1(C), higher values of ANG indicates the larger extent of departure of α9 from the BCG. Along this axis, the most notable characteristic of the plot for ABX [Figure 2(A)] is the presence of a steeper barrier (~ 6-8 kcal/mol) around BCG boundary (ANG ~ 95 - 105˚) that separates two valleys (local and global minima), whereas rise of the 10 ACS Paragon Plus Environment

Page 10 of 40

Page 11 of 40 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

Journal of Chemical Information and Modeling

landscape height is much gentle for BBX [Figure 2(B)] while α9 moving away from the pocket. There lies a relatively larger barrier (~ 8-10 kcal/mol) outside the BCG (ANG ~ 100˚ - 110˚). Beyond the barrier, though both landscapes become nearly horizontal along ANG axis, corresponding states of ABX and BBX markedly differ in energy. This will be more comprehensible from the discussion of 1D projections shown in Figure 3. Comparison of the landscapes within BCG (ANG ~ 80˚ - 100˚) shows that slopes of the surfaces along the E2ED axis on the two sides of the valley are relatively steeper in the case of BBX. This indicates the narrower distribution of E2ED (~ 26 Å – 28 Å) inside the groove (ANG ~ 80˚ 100˚) correspond to lower energy conformers (indicated by blue shed) in the case of BBX compared to ABX (blue shed encompasses E2ED ~ 26 Å – 30 Å). In contrast to flexible (E2ED ~ 26 Å – 30 Å) α9 conformers of ABX at the global minima (inside the BCG), a broad-flat local minimum (green shade) lies outside (i.e. in bulk aqueous phase) that correspond to only compact structures (E2ED ~ 26 Å – 28 Å), but no such comparable local minima exist for BBX. In fact, the overall surface is much flatter and lower in energy outside the groove for ABX in comparison to a rougher surface of BBX. In Figure 3, the 2D PMF was projected along ANG at fixed values of E2ED (using Eq. 2 described under Methods section). This was done to locate the barrier position in reference to BCG boundary (ANG ~ 100˚) and to get the idea of possible routes of α9 exclusion in terms of defined RCs. Figure 3(A) shows ABX landscapes at varying E2ED differing only in their barrier height but identical in barrier position (ANG ~ 100˚) and shape of the curve. Moreover, comparable ground states of activation (α9 bound states; ANG ~ 80˚) at different E2ED separated mostly by 2-4 kcal/mol (solid lines) suggest presence of a flexible C-tail (α9) bound to the BCG with a wider range of helicity [Figure S4(A)] in ABX. Helicity was estimated using

11 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

‘alpha’ collective variable49 implemented in NAMD. The ‘alpha’ collective variable deals with an angle scoring function considering backbone Cα atoms and a scoring function for the backbone H-bonds (SI text). For BBX [Figure 3(B)], ground states comparable with ABX correspond to only compact helical conformations (E2ED ~ 26 – 28 Å) (solid lines). In fact, BBX α9 maintains its helicity even at higher E2ED [Figure S4(B)] within the BCG that suggests presence of an ‘extended helical’ conformation instead of a ‘lower helical’ termini like ABX at higher E2ED. As it is mentioned earlier that ABX conformations with various length of α9 are of comparable energies (differing by 2-4 kcal/mol), those are expected to be of comparable populations too. In consequence, BIM can notably bind to any conformation representing those populations with similar probability. In Figure 4, projections are compared (1:1) for ABX and BBX at different E2ED and plausible routes (energetically favorable route on top of energy surface) are suggested for α9 exclusion (solid cyan lines) upon BIM binding to apo conformations. Unanimously all the projected PMF profiles of ABX and BBX appears to be intersecting around BCG boundary [vertical black dashed lines in Figure 4(A to F)] which means at this point systems can choose its course of action regarding removal of α9. Beyond this point, ‘open’ states (i.e. α9 excluded from the BCG) corresponding to ABX are of lower free energy (for both compact and extended helices) in comparison to those of BBX. Hence all the suggested pathways (SP: cyan line) indicated that once α9 crosses the BCG boundary, profile of the ABX would become a thermodynamically preferred track and therefore BIM may get disengaged immediately from BBX to follow the low energy track. Thus, comparison of PMF profiles provides a thermodynamic insight into the experimentally suggested ‘hit-and-run’ model where transient encounter of activator BH3-only proteins (e.g. tBID, BIM etc.) was proposed to be sufficient for pro-apoptotic BAK and BAX activation rather than formation of a stable pro-

12 ACS Paragon Plus Environment

Page 12 of 40

Page 13 of 40 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

Journal of Chemical Information and Modeling

apoptotic: activator complex.10, 19, 24, 66 In addition, as BIM unbiasedly destabilizes the ground states (ANG ~ 80˚) (whether compact or extended α9), it reduces the free energy difference between the α9 bound and released states considering ‘transient’ binding of BIM (i.e. energy difference between global (inside the groove) and local (outside the groove) minima along SP in Figure 4). This reduction is ~ 1-2 kcal/mol for compact α9 [Figure 4(A) – (C)] but ~2-4 kcal/mol for extended helices [Figure 4(D) – (F)]. Such favor for extended helices is additionally benefited from reduction in barrier heights. ABX and BBX show comparable barriers (∆∆E ~ 0.5 – 1 kcal/mol) for compact α9 conformations (E2ED ~ 25 – 28 Å), whereas BBX yields smaller barriers (lowering by ~3–4 kcal/mol) for extended helical conformations (E2ED > 28Å) in comparison to that of ABX. So, the population of extended helical α9 seems to be triggered more favorably than that of compact helical α9. Calculation of Lowest Free Energy Path (LFEP)

67

provides transition routes only between two minima in the 2D energy landscape (Figure S3). That’s why in the present study, LFEP overlooks the possibility that BIM can bind to a flexible range of ABX population with similar probability. Later 1D projections along LFEP [Figure 3(C)] also confirm that α9 exclusion is energetically comparable (~8 – 8.5 kcal/mol) among ABX and BBX for compact helical α9 [similar to Figure 4(A)-(C)]. So, if activation follows LFEP, the only difference in the landscape is the position of the barrier which has clearly been shifted outside the groove for BBX in comparison to ABX that may kinetically favor α9 of BBX over that of ABX to cross the BCG boundary.

Contact analysis: The MDABX and MDBBX trajectories (200ns each) were analyzed to sort out important contacts that changed their life-time remarkably from one system to another.68 To simplify the complex configurational space, the strategy by Johnson et al.69 was adopted. Here, a distance

13 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

cutoff of 4.2 Å between any two heavy atoms of two residues during simulation is considered as accountable ‘contacts’.70,

71

Later, the difference in the probability of a contact formation is

calculated as ∆f = (fBBX– fABX) where f is defined as followsf =

Number of frames the contact lasts Total no of frames

Though f can vary from 0 (contacts that never formed) to 1 (contacts intact throughout the simulation), contacts with transition stability [0.1 < f < 0.9] are keys to contribute differences in conformational ensembles.69,

70

Parameter ∆f sorts out the contacts differing in their lifetime

between the two systems (i.e. ABX and BBX), where a positive value (∆f > 0) for a contact pair implies that the contact is more stable in BBX and a negative value (∆f < 0) implies that the contact is more stable in ABX. Figure 5 elucidate those contacts which are either dominant [Figure 5(A)] or unique [Figure 5(B)] in any of the two systems. As only distance criteria are used to build this map, it is somewhat reluctant in regard to ‘nature of the contacts’ (i.e. electrostatic or hydrophobic). Therefore, only a few of the clusters in the map can be argued to be important in terms of specific ‘nature of the contact’ or agreement with previously available experimental data can also be checked. From Figure 5(A), α1-α5 contacts [region-(i)] appear to be ~60% (i.e. average ∆f ~ 0.6) more stable in BBX trajectory, whereas α5-α6 contacts [region-(iv)] are ~30% more stable in ABX (i.e. average ∆f ~ -0.3). These contacts are graphically shown in Figure S5. This difference can be attributed to the binding of BIM in BBX and it implies that BIM binding is enhancing the stability of the α1-α5 contacts (core) and loosening the α5-α6 contacts (core-latch). This is in agreement with the “core-latch separation” hypothesis for membrane-embedded activated monomers proposed by Czabotar et al.

26

. Their proposal in regard of triggering the core-latch

separation by loosening of packing between Ile19 of α1 and Ile149, Ile152 of α7 is also 14 ACS Paragon Plus Environment

Page 14 of 40

Page 15 of 40 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

Journal of Chemical Information and Modeling

reconciled from simulation data. Those proposed contacts are found to be ~20% more stable in ABX compared to BBX (i.e. average ∆f ~ -0.2). A higher population of relatively stable ABX contacts (green dots) in region-(ii) (Figure 5A) suggests less stable contacts between α9 and BCG residues upon rear pocket BIM binding in the BBX. This suggests that α9 is prone to get exposed to the solvent in case of BBX, at least more than the same for ABX, which is the primary requirement of ‘activation’. Experimental report72 suggested that Leu70 (α2) is a key BCG residue that interacts with α9 residues, e.g. Thr174, Val177, and Ala178 in ABX. In agreement with that, Leu70 – Ala178 and Leu70 – Leu181 interactions are found to be ~15% and ~30% more stable in ABX in comparison to the same in BBX (i.e. ∆f ~ -0.15 and -0.3 respectively) (Figure S6). Such reconciliation of experimental findings ensures that conformations from ABX and BBX ensembles are a considerable representation of ‘inactive’ apo BAX and ‘activated’ BIM bound BAX respectively. There are dominant BBX contacts that involve α3-α4 in one part and α5-α7 in another part [region-(iii)]. Relatively long-lived hydrophobic contacts, such as Val91 (α4) - Leu120 (α5) (~30% more stable; ∆f ~ 0.3), Thr85 (α3-α4 loop) - Ala124, Val129 (α5-α6 loop residues) (~20% more stable; ∆f ~ 0.2) suggest stable core contacts in the complex. Moreover, there were several interactions around BCG, involving an aromatic ring and charged side chains (i. e. potential cation-π interacting pair) those have changed significantly upon apo to complex transition, e.g. Phe105 (α4-α5) – Arg147 (α6-α7) get notable stability (~70%; ∆f ~ 0.7) in BBX. Although the cation-π interactions are not explicitly modeled in molecular mechanical force fields, the stabilizing effect is implicitly included in the parameters to some extent. In the force field, for a phenyl ring, the negative charged on the carbon and equal positive charges on carbon results in quadruple moments and therefore when a positive charge (e.g. Arg, Lys side chains) approaches

15 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

towards the ring it can enjoy some stability until the van der Walls repulsion becomes active. There are reports available in the literature that underscores the roles of such interactions in stabilizing the folded structures of proteins. For example, previous study has revealed that there is a substantial number of structures available in Protein Data Bank where Arg and Trp are biased in cation-π interactions.73 Generally, a 6Å distance cut-off is used to identify cation-π interaction pair. Likewise, as indicated by the change in the most probable distance between Trp170 (α9) side-chain ring and Arg65 side chain (α2) [Figure 5(D)], these two residues appear to be a potential cation-π interacting pair in the ABX trajectory. On the other hand, BIM binding altered such interaction by recruiting Arg65 into a strong salt-bridge interaction network involving polar α2 (Arg65, Lys64, Asp68 and Asp71) and α1 residues (Asp33, Arg34, and Arg37) [Figure 5(C)] and eventually that cation-pi interaction was lost in the BBX trajectory. These polar contacts are reflected as region-(v) of Figure 5(A). The pairwise interaction energy of participating residues in this salt-bridge network was averaged over the trajectory and plotted in Figure S7. Though in ABX, Asp33, and Arg34 of α1 are participating in a few specific pair of interactions (e.g. Asp33 - Lys64, Arg34 – Asp68), Arg37 is completely mute in response to any interacting partner because it is in a loop (see Methods section). On the other hand, Arg37 is the part of a helix in BBX, which makes all three residues (Asp33, Arg34, and Arg37) involved in substantial interaction with α2 residues. The newly formed contacts involving Arg37 are also evident from region-(v) of Figure 5(B) where exclusive contacts with substantial stability (∆f ~ 0.4) are shown. Network of most dominant contacts (∆f > 0.4) discussed in Figure 5(A) is presented in Figure S8(A, B). Now it is comprehensible from the diagram that the number of contacts involving α1, α5, and α6 has been reduced significantly during complex formation.

16 ACS Paragon Plus Environment

Page 16 of 40

Page 17 of 40 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

Journal of Chemical Information and Modeling

Figure 5(B) shows the newly formed as well as completely lost contacts in BBX trajectory. Most notable new contacts (∆f > 0.73) are hydrophobic in nature [region- (iii) and (iv)], e.g. Met38 (α1) with (Leu122, Cys126) of α5 and Trp107 (α5) with Trp151 (α7) etc. However, several significant polar contacts (∆f ~ -0.8) involving α1 are lost upon BIM binding [region-(i) and (ii)], e.g. Ser16, Glu17 (α1) – Arg145 (α6) ; Arg34 (α1) – Tyr115 (α5) ; Arg34 (α1) – Lys119 (α5) etc. Moreover, a substantial number of contacts involving α1–α2 loop and rear pocket residues have been lost in the transition of ‘close’ to ‘open’ loop conformation (Figure S9) upon BIM binding [region- (vi) of Figure 5(B)] in agreement with Gavathiotis et al.19 and our initial modeling (see Methods). It offers a possibility that the flexible ‘open’ conformation of the α1–α2 loop may have an effect on the dynamics of α9 outside the BCG and consequently, that would influence further activation steps which need a separate study to explore. Similar to Figure S8(A, B), network of unique contacts is shown in Figure S8(C, D). Thus BIM binding to the rear pocket exposes α9 as well as provoke ‘core-latch separation’ in cytosolic phase. It also affects BC groove dynamics by complimenting between cation-π and salt-bridge interactions.

Principal Component Analysis: How BIM induced rearrangement of the interaction among BAX residues affects its global dynamics, can be identified by extracting the principal modes of motion from the equilibrium trajectories of ABX and BBX using principal component analysis (PCA). PCA was performed on the MDABX and MDBBX trajectories, which revealed contribution of the first 3 PC modes (PC1-PC3) for ABX to be 26.35%, 9.02%, and 6.51% to the total variance and the same for BBX to be 30.57%, 9.78%, and 6.64% respectively i.e. these top PCs cover almost 50% of the overall 17 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

motion. For comparing the conformational spaces sampled for ABX and BBX, another set of PCs were calculated on a single comprised dataset of both ABX and BBX frames. Projection of conformations from each of the systems (MDABX and MDBBX) onto two top PCs shows a narrower distribution of ABX [Figure 6(A)] ensemble in comparison to BBX [Figure 6(B)]. This difference in sampling space is also evident from the calculation of cosine correlation between first 3 PC modes of ABX and BBX systems (Table-S3). Observed overlaps are ~30% among PC1 modes and ~46% among PC2 modes which individually contributed ~30% and ~10% to the global variance respectively. This infers that BIM significantly perturbed the global dynamics (in terms of top PCs) of BAX monomer. It is important to mention that all the PCs are calculated for the BAX monomer from both ABX and BBX trajectories. Looking for substantial structural modifications upon BIM binding is important to probe for the difference observed in the sampling space. Square fluctuation of residues along PC1 [Figure 6(C)] has been calculated as that contributes maximum to the displacement along the trajectory. It shows a much higher fluctuation around N-terminal of α6 (residues 130 to 140) for BBX trajectory compared to ABX because α6 has lost its helicity significantly upon interaction with BIM at the rear pocket (made up with α1 and α6). Later, conformations were projected on a helicity (α6) vs. PC1 surface [Figure 6(D) and (E)] and time evolution is shown on a color scale of red (starting of trajectory) to blue (end of trajectory). In comparison to BBX, ABX [Figure 6(D)] ensemble seems to be populated by a lower range of PC1 values (-1.5 - 0.5) with helicity greater than ~70% (i.e. helicity is mostly maintained throughout the trajectory). Whereas for last 50ns, BBX conformers [Figure 6(E)] appeared in the range of lower helicity (less than 60%) and higher PC1 values (-0.5 – 1.0). Apart from α6, residues of the C-terminal tail (α9) are also relatively more fluctuating in BBX system [Figure 6(C)]. Further, BC groove volume has been

18 ACS Paragon Plus Environment

Page 18 of 40

Page 19 of 40 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

Journal of Chemical Information and Modeling

calculated using POVME_2_0_1 package74 and compared [Figure 6(F)]. The average pocket volume for ABX is larger (534 ± 37 Å3) than BBX (460 ± 43 Å3). One point worth mentioning is that the magnitude of the calculated volume may vary from one algorithm to another due to differences in the definition of a pocket and in the probing parameters. Therefore the usage of pocket volumes for qualitative comparison is more meaningful than banking on their exact magnitudes.

Network analysis on Equilibrium Trajectories: To unravel the information (in terms of interaction) transfer pathways among various domains of BAX, time-averaged interactions between residues were used to elucidate networks and to group them into different communities (see Methods). Important residues belong to each community have been listed in Table-S4. ABX shows 9 communities, referred as C-1, C-2 etc. up to C-9, with 1 isolated community (C-2) containing α1-α2 loop residues, whereas 10 communities (C-1 to C-10) were found for BBX with 2 isolated communities (C-1 and C-3) comprised of only 3 residues (Figure 7). In the case of ABX [Figure 7(A)], α9, as well as BCG (i.e. α2 to α5) residues, appear to be distributed over 3 specific communities (C-7, C-8, and C-9). As intra-community nodes interact strongly, dynamics of α9 can be argued to be strongly correlated with dynamics of BCG in ABX system. Inter-community connectivity between a community with α9 dominance [C-9: 70% of α9 residues] and communities with major BC groove residues (C-7, C-8, and C-4) also suggests stronger correlation [Figure 7(B)]. In these community maps [Figure 7(B) and (D)], the width of an edge connecting two communities is proportional to the sum of betweenness (number of the pairwise optimal path that crosses that edge) of edges connecting nodes of those two

19 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

Page 20 of 40

communities. The major dissimilarity of BBX with apo system lies in intra-community membership distribution. As opposed to ABX, here 90% residues of α9 belong to a single community (C-7) along with few α2-α3 residues [Figure 7(C)]. Yet C-7 is connected to C-8, C10 and C-6 [Figure 7(D)] where C-6 and C-8 contain few BCG residues. So BIM primarily isolates residues of α9 from colonies of BCG residues and redistributes them into a single community. Still, the question remains – how does BIM interact with α9 at the BC groove upon transient binding at the rear groove? Precisely, if they cannot interact directly, they can cross-talk. To dissect that, sets of residues from BIM and α9 were chosen first and then shortest (optimal) paths were calculated (details are described in Methods) among all possible pairs of those residues (Figure 8). From BIM, residues Arg153, Ile155, and Asp157 were chosen as experiment

25

suggested those to be critical for BIM bound BAX activation. From α9, residues facing the hydrophobic BC groove and crucial for α9 – BCG interaction, were considered (see Table S5). It was found that paths connecting BIM with Thr174, Val177, and Ala178 of α9 are mostly conserved through residues from α1 and α2 [Figure 8(A)]. Most frequently appearing route is: Thr174 (α9) – Ile66 (α2) - Lys64 (α2) - Phe30 (α1) - Gln28 (α1) - Phe159 (BIM), whereas Leu70 (α2) - Gly67 (α2) - Phe30 (α1) - Gln28 (α1) - Phe159 (BIM) is second in the rank (Table S5). The sequence Phe30 (α1) - Gln28 (α1) - Phe159 (BIM) is common in all the routes. Other than this α1-α2 mediated path, there exists an alternative set of conserved residues belong to α6 and α4 (reverse face of α1 and α2) which connect BIM with residues near C–terminal of α9 (Leu181, Thr186, and Lys189) [Figure 8(B)]. Lys189 and Leu181 were reported to be important for BIM bound activation.19 Here most consistently occurring sequence is Leu185 / Trp188 (α9) – Val95 / Arg94 (α4) – Phe93 / Phe92 (α4) – Trp139 (α6) – Leu141 (α6) – Leu152 (BIM) (see Table S5). 20 ACS Paragon Plus Environment

Page 21 of 40 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

Journal of Chemical Information and Modeling

Among residues in this sequence, contacts like Leu185 (α9) – Val95 (α4) and Leu141 (α6) – Leu152 (BIM), are also observed to be contacts with highest edge connectivity for C-7:C-6 and C-6:C-4 respectively (Figure S10) where C-4 of BBX contains the whole BIM peptide, C-7 contains the majority of α9 and one of the cross-talk routes is via C-6.

Discussion: Though there is convincing evidence7,

19, 25

in favor of rear pocket BIM triggered BAX

activation, how this pivotal step initiates the process by displacing α9 from BC groove in the first place is vaguely understood.26 Study of the dynamics and energetics certainly provides assistance to bring some clarity in the associated mechanism. The 2D energy surfaces (Figure 2) clearly give the impression that BIM binding at the rear pocket modified the landscape of α9 release from the BCG. Not only the slope near the BCG boundary is changed but also the surface becomes rougher and energetically higher beyond the groove upon BIM binding. It suggests that states corresponding to α9 excluded from the BC groove are more favorable for ABX in comparison to the BBX system. Later, comparison of 1D PMF projections (Figure 4) suggests pathways where BIM would prefer to leave the BBX system once α9 leaves its native residence (i.e. the BCG), switching the system into the track of ABX. This explains the thermodynamic requirement for ‘transient binding’ of activator protein which is sufficient to activate BCL-2 family of proteins.10,

19, 24, 66

Moreover, as BIM unanimously reduces free energy difference

between the local and global minima along SP (Figure 4) for compact or extended helices, it implies that whatever may be the conformation of α9 within the BCG, the decision of its release 21 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

Page 22 of 40

is solely determined by the ‘transient’ mode of activator (i.e. BIM) binding. In addition, such activation is more pronounced for the extended helical α9 conformers due to an additional benefit from lowering of barrier heights in comparison to compact helical α9 [Figure 4(D) – (F)]. As landscapes within the BCG (i.e. ANG < 100˚; Figure 2 and Figure 3) suggest more compact α9 conformations (E2ED~ 26 Å – 28 Å ) inside BCG of BBX compared to that of ABX (E2ED ~ 26 Å – 30 Å), it is likely that α9 has to maintain a compact shape to remain bound inside the BCG of BBX, which is in agreement with previous reports that claimed α9 to be more flexible in apo BAX compared to the same in the BAX:BIM complex.17,

39

To this end, it can be

qualitatively argued that reduced number of microstates (in terms of flexible α9 conformers) implies a relatively lower probability of finding the α9 inside BCG for BBX system compared to ABX. In a nutshell, BH3 encounter (BIM in this case) at the rear side of the BC groove has some allosteric control over the compactness of α9 packed in the groove. The next step would be to look at the residue level contacts in BAX monomer which get modified upon BIM binding (Figure 5). In line with observations from energetics, Contact analysis on equilibrium trajectories (i.e. MDABX and MDBBX) shows weaker contacts of α9 with BC groove residues in BBX system [region-(ii) of Figure 5(A)]. This is also reflected in Figure 6(C) where residues of α9 are more fluctuating along PC1 in BBX compared to ABX, but RMSF (Figure S2) of the trajectories doesn’t show any such differential fluctuation in α9 residues. This suggests that weaker contacts between α9 and BCG residues have a notable impact on the most contributing mode of the global dynamics (PC1). BIM binding also disrupted the cation-π interaction between Trp170 (α9) and Arg65 (α2), by indulging Arg65 into a salt-bridge interaction network involving polar residues (Asp33, Arg34, and Arg37) of α1 and α2 (Arg65, Lys64, Asp68, and Asp71). In a previous computational study,39 charged residues of the BH3 22 ACS Paragon Plus Environment

Page 23 of 40 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

Journal of Chemical Information and Modeling

domain (α2) i.e. Lys64, Asp68 and Asp71 were reported to maintain interactions with α1 upon BIM binding. Mutagenesis studies also confirmed their role in dimerization upon activation.66 Moreover, Asp68 and Asp71 are reported to be important for maintaining the soluble monomeric state of BAX.72 Therefore, disruption of the native cation-π interactions involving α9-BCG residues and subsequent formation of a salt-bridge network between α1 and α2 upon BIM binding is arguably critical for BAX activation. Due to the interaction with BIM at the rear pocket, α6 loses its helicity [Figure 6(C, E)] quite significantly around N-terminal and that causes loosening of α5-α6 contacts (core-latch contacts) [region-(iv) of Figure 5(A)]. Yet several new α1-α5 hydrophobic contacts (∆f > 0.73) are formed [region- (iii) and (iv) of Figure 5(B)] and few get significant stability (∆f ~ 0.60) [region-(i) of Figure 5(A)] upon BIM binding which suggest for stable core contacts. These observations emphasize the ‘core-latch separation’ hypothesis proposed by Czabotar et al.26 where dislodgement of the core (α1-α5) from the latch (α6-α8) domain was suggested to be an important structural change for a membrane-embedded activated monomer, which supposedly processed by decoupling of anti-parallel α5-α6 helices. Taken together experimental and these in-silico observations, it can be argued that core-latch separation initiates the primary stage of structural modification after BIM binding in the cytosolic phase and that further gets matured upon membrane integration. So experimentally observed different steps of activation,7,

19, 25, 26

are found to be concurrent with this study of

dynamics. PCA was performed on equilibrium trajectories to understand the essential dynamics of the systems. Relatively broader distribution of BBX on top PCs indicates greater exploration (in terms of more flexible conformers) of the conformational space upon BIM binding [Figure 6(A, B)]. Significant differences observed from contact map analysis (discussed earlier) clearly hint 23 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

for distinct conformational ensembles of BAX monomer for ABX and BBX systems and that is what reflected in top PC modes. Later, comparison of conformations along PC1 shows that loss of α6 helicity causes higher values of PC1 for BBX i.e. more positional fluctuation along this PC mode [Figure 6E]. So, α6 can be crucial for causing a difference in sampling space between ABX and BBX trajectories. Structural effect of BIM binding on the BC groove is directly evident from reduced BC groove volume upon BIM binding [Figure 6(F)]. Network Analysis was performed to get the idea of how different regions of the BAX monomer communicate with each other and how that connectivity is modified upon BIM binding. Community analysis reveals a stronger correlation between α9 and BCG residues in ABX compared to BBX [Figure 7(A, B)]. Due to lack of intra-community membership, BIM primarily weakens communication of α9 with BCG which gets compensated by inter-community interactions to some extent. As intra-community nodes interact strongly than inter-community nodes, this holds the remark that BIM exposes α9 due to weaker interaction with BC groove. In the concern of how does BIM exert its influence on α9, shortest paths were calculated among residues of α9 and BIM [Figure 8]. Two probable conserved allosteric pathways from BIM are explored, where one branches through α1 and α2 to connect residues near the N-terminal end of α9 (i.e. Thr174, Val177 and Ala178) and another through α6 and α4, critical for C–terminal end residues (i.e. Leu181, Thr186, and Lys189). It seems that α9 can be accessed by BIM from both the faces (α1-α2 and α6-α4) due to the symmetric globular shape of the protein core. Previous computational study 39 also proposed allosteric communication through residues of α6 and α4 yet specific residues reported are not in complete agreement with the present study. Suggested critical residues [Gln28 (α2), Gly67 (α2), and Thr174 (α9)] for BIM induced BAX activation from NMR experiments19,

25

are in complete agreement with our calculated paths. Other α2 24 ACS Paragon Plus Environment

Page 24 of 40

Page 25 of 40 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

Journal of Chemical Information and Modeling

residues appeared in the path, Leu70 and Ile66, alongside α9 residues Val177 and Ala178, were discussed important for the subcellular distribution of BAX monomer.72

Conclusion: The molecular mechanism of BAX activation through rear pocket triggering has been presented in this report in the light of its activation energy profile. All-atom MD simulation has probed and revealed how BIM kicks out α9 from a distal site by recalibrating the associated thermodynamic profile. The investigation has also provided a benchmark for the choice of E2ED and ANG as two possible reaction coordinates for expressing PMF of such process, particularly for the similar type of systems. Success of the choice has provided confidence for further investigation to understand the complex molecular mechanism of activation process related to other proteins of BCL-2 family. The PMF landscapes have shown that BH3 encounter (BIM in this case) at the rear side of the ligand binding BC groove controls the flexibility of α9 packed within the groove and modulates the landscape of α9 exclusion depending on ‘which’ population (compact or extended helical α9) of apo BAX is getting triggered. The systematic comparison of 1D PMF for apo and complex systems provides the thermodynamic basis of the well-established ‘hit-and-run’ model of activation which reconciles the previously reported hypothesis from experimental data.10, 19, 24, 66 Analysis of the contact map has shown significant alteration in a network of saltbridge and hydrophobic clusters of interactions which are crucial for α9 exposure upon BIM binding. Two potential allosteric routes have been revealed from this investigation which suggests that the induction effect of BIM can propagate to BC groove through α1-α2 as well as through α6-α4 helices i.e. both faces of the protein are available for BIM to allosterically access α9 from the rear pocket. So far the authors are aware, for the first time, the effect of rear pocket BIM binding on the activation energy landscape of BAX is hereby delineated. 25 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

26 ACS Paragon Plus Environment

Page 26 of 40

Page 27 of 40 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

Journal of Chemical Information and Modeling

Figures:

Figure 1: Cartoon representation of BAX monomer structure: (A) Apo BAX (ABX) structure. Helix-1(i.e. α1) to Helix-9 (i.e. α9) is shown in rainbow color scale. (B) BAX:BIM complex (BBX) structure where BAX monomer is shown in surface representation. BIM peptide (Gray) and α9 (Red) helix are represented in cartoon structure. Rear pocket or BIM interaction site is shown with the cyan surface whereas α9 binding BC groove is shown with the light yellow surface. Description of Reaction coordinates (RC): (C) Angle RC (alpha carbon of Ile66, Thr172 and Gly179 are referred as i, j and k respectively). (D) End-to-end distance (backbone carbon of Trp170 and Trp188 are referred as p and q respectively). 27 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

Figure 2: Free energy profile of α9 release from BC groove (BCG) in (A) apo BAX (ABX) and in (B) BAX:BIM complex (BBX). Energy values are mapped with color gradient. Blue to red shed describes regime of lower energy to higher energy. Regions of the barriers are pointed by red arrows. Representative structures of specific regions are provided and α9 is colored in light green. The distance between orange spheres is representing end-to-end distance (E2ED) RC.

28 ACS Paragon Plus Environment

Page 28 of 40

Page 29 of 40 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

Journal of Chemical Information and Modeling

Figure 3: 1D Projections of 2D PMF along angle RC (ANG) at different E2ED for (A) ABX and (B) BBX. Profiles having ground state (ANG ~ 80˚) energy in the range 6 - 10 kcal/mol and greater than 10 kcal/mol are represented by solid lines and dotted lines respectively. (C) 1D projections along LFEP. The vertical dashed lines around ~100˚ represent boundary of the BC groove.

29 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

Figure 4: (A-F) represent 1:1 comparison of projections between ABX (solid lines) and BBX (dotted lines) systems at different E2ED. Cyan solid lines represent probable suggested pathways for α9 removal upon BIM binding at the rear pocket. Colors representing projections at each E2ED is chosen from Fig-3A. The black vertical dashed lines around ~100˚ represent boundary of the BC groove (BCG). The global minima inside the BC groove (ANG < 100˚) and local minima outside (ANG > 100˚) along the suggested pathway are marked by black rectangle and position of energy barrier along the suggested pathway is marked with a black dot. Barrier

heights along the ABX and BBX pathways are provided for each projection.

Figure 5: Mapping of dynamically significant contacts those have changed upon BIM binding at the rear pocket. (A) Differences in the fraction of time a contact lasts in any of the two trajectories (ABX and BBX) are mapped. Dots with green shades (-ve ∆f values) represent contacts more stable in ABX trajectory and red shaded dots (+ve ∆f values) represent stable contacts in BBX. (B) Contacts unique to a trajectory are mapped. Greenish dots are representing contacts that last at least 10% of the total simulation time in ABX trajectory but was completely 30 ACS Paragon Plus Environment

Page 30 of 40

Page 31 of 40 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

Journal of Chemical Information and Modeling

lost in BBX trajectory. Similarly, reddish dots represents contacts newly formed in BBX (lifetime > 10%) but was missing in ABX. (C) Representation of the observed salt bridge network involving α1 and α2 polar residues for ABX (left panel) and BBX (right panel) systems respectively. (D) Distribution of the separation between Trp170 and Arg65. Blue colored curve is for ABX system and the red curve is for BBX system. ABX has fairly sampled a favorable cation-π distance with the peak at ~5.2 Å, whereas BBX distribution peaked at ~9Å which is far beyond the favorable cation-π distance (~6 Å).

Figure 6: Projection of (A) ABX and (B) BBX conformations onto 2 top PCs. Time evaluation is projected on a black to yellow color scale. (C) Square fluctuation of residues along PC1. The flanking, disordered loops (residues 10 – 17 and 36 - 54) are not considered to avoid larger scaling in y-axis. The inset diagram represents a superposition of structures of the last frame of respective equilibrium trajectories where orange circle indicates the region of α6 that lost its helicity in BBX. The same region is also marked in the RMSF plot with the orange circle where BBX is showing relatively larger fluctuation. (D, E) Projection of ABX and BBX trajectories on a Helicity (α6) vs PC1 surface. Time evolution of frames is shown in red to blue color scale. (F) Time evolution of BC groove volume for ABX and BBX trajectories.

31 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

Figure 7: Community analysis of the network based on equilibrium simulation trajectories. (A) ABX and (C) BBX systems are colored according to calculated community membership. Community colors along helices of the monomer are presented with color bars. (B, D) Community network representation for ABX and BBX respectively. Communities are colored according to (A) and (C). Communications between major α9 containing community (blue block) and other communities are shown using blue arrows. Width of the edges is proportional to

32 ACS Paragon Plus Environment

Page 32 of 40

Page 33 of 40 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

Journal of Chemical Information and Modeling

the strength of interactions between connected communities. Isolated communities are not shown in community networks.

Figure 8: Shortest paths of communication between BIM (at the rear pocket) and α9 (at the BC groove) in BBX. Orange spheres represent residues considered for shortest path calculation and silver spheres represents residues through which the shortest paths connect BIM and α9 residues i.e. orange spheres. Sticks of different colors represent edges connecting two adjacent nodes (Cα of residues) in a path. (A) Conserved paths through α1 and α2 residues that connect BIM with Ala178, Val177 and Thr174 of α9. (B) Conserved paths through α6 and α4 residues that connect BIM with Leu181, Thr186 and Lys189 of α9. In both the routes, same set of BIM residues i.e. Arg153, Ile155 and Asp157 (orange spheres in BIM) are considered as terminal residues of the shortest path calculation.

33 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

ASSOCIATED CONTENT Supporting Information: The Supporting Information is available free of charge. Table-S1: Range of residues in each helix of BAX monomer; Table-S2: List of trajectories; Figure S1 – S2: Equilibrium simulation stability check; Figure S3: LFEP calculation; Figure S4: Helicity of α9 projected on 2 chosen reaction coordinates; SI Text: Description of alpha collective variable; Figure S5: Graphical representation of α1-α5 and α5-α6 hydrophobic packing; Figure S6: Position of Leu70 – Ala178 and Leu70 – Leu181 interactions in the differential contact map and subsequent graphical representation; Figure S7: Pairwise average interaction energy among polar residues of α1 and α2; Figure S8: Network view of dynamically significant contacts those have changed upon BIM binding; Figure S9: Transition of ‘closed’ to ‘open’ conformation of the α1-α2 loop; Table-S3: Cosine correlation of PC modes; Table-S4: Residue distribution among communities of ABX and BBX networks; Table-S5: Shortest paths among BAX α9 residues and BIM residues; Figure S10: Cross-talk pathways between communities; SI Text: Validation of convergence of PMF profiles; Table-S6: Standard deviation of system force for checking of convergence of PMF; Figure S11: Distribution of ANG RC from short simulations performed for committor analysis; Figure S12: Distribution of committor (PDF).

This information is available free of charge via the Internet at http://pubs.acs.org

AUTHOR INFORMATION Corresponding Author

34 ACS Paragon Plus Environment

Page 34 of 40

Page 35 of 40 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

Journal of Chemical Information and Modeling

*Email: [email protected]

ORCID Souvik Sinha: 0000-0003-4964-7713 Atanu Maity: 0000-0002-7795-692X Shubhra Ghosh Dastidar: 0000-0002-0692-8139 Notes The authors declare no competing financial interest.

ACKNOWLEDGEMENT SS thanks UGC (University Grant Commission), Government of India for funding. AM thanks CSIR (Council of Scientific and Industrial Research), Government of India for funding. Authors acknowledge Dr. M Moradi (U. Arkansas) for the LFEP code. ABBREVIATIONS Cyt C : Cytochrome C ; MOM : Mitochondrial outer membrane ; BCG : BC-groove ; MD : Molecular dynamics ; ABX : apo BAX ; BBX : BAX:BIM complex ; PMF : Potential of Mean Force; ABF : Adaptive Biasing Force ; RC : Reaction Coordinates ; ANG : Angle of opening (Ѳ) reaction coordinate; E2ED : End-to-end distance ; PCA : Principle component analysis ; LFEP : Lowest Free Energy Path

References:

35 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

1. Danial, N. N.; Korsmeyer, S. J. Cell death: Critical control points. Cell 2004, 116, 205-219. 2. Leber, B.; Lin, J.; Andrews, D. W. Embedded together: the life and death consequences of interaction of the Bcl-2 family with membranes. Apoptosis 2007, 12, 897-911. 3. Youle, R. J.; Strasser, A. The BCL-2 protein family: opposing activities that mediate cell death. Nat. Rev. Mol. Cell Biol. 2008, 9, 47-59. 4. Cheng, E. H.; Wei, M. C.; Weiler, S.; Flavell, R. A.; Mak, T. W.; Lindsten, T.; Korsmeyer, S. J. BCL-2, BCL-X(L) sequester BH3 domain-only molecules preventing BAX- and BAK-mediated mitochondrial apoptosis. Mol. Cell 2001, 8, 705-711. 5. Wei, M. C.; Zong, W. X.; Cheng, E. H.; Lindsten, T.; Panoutsakopoulou, V.; Ross, A. J.; Roth, K. A.; MacGregor, G. R.; Thompson, C. B.; Korsmeyer, S. J. Proapoptotic BAX and BAK: a requisite gateway to mitochondrial dysfunction and death. Science 2001, 292, 727-730. 6. Cartron, P. F.; Oliver, L.; Mayat, E.; Meflah, K.; Vallette, F. M. Impact of pH on Bax alpha conformation, oligomerisation and mitochondrial integration. FEBS Lett. 2004, 578, 41-46. 7. Kim, H.; Tu, H. C.; Ren, D.; Takeuchi, O.; Jeffers, J. R.; Zambetti, G. P.; Hsieh, J. J.; Cheng, E. H. Stepwise activation of BAX and BAK by tBID, BIM, and PUMA initiates mitochondrial apoptosis. Mol. Cell 2009, 36, 487-499. 8. Kim, H.; Rafiuddin-Shah, M.; Tu, H. C.; Jeffers, J. R.; Zambetti, G. P.; Hsieh, J. J.; Cheng, E. H. Hierarchical regulation of mitochondrion-dependent apoptosis by BCL-2 subfamilies. Nat. Cell Biol. 2006, 8, 1348-1358. 9. Kuwana, T.; Mackey, M. R.; Perkins, G.; Ellisman, M. H.; Latterich, M.; Schneiter, R.; Green, D. R.; Newmeyer, D. D. Bid, Bax, and lipids cooperate to form supramolecular openings in the outer mitochondrial membrane. Cell 2002, 111, 331-342. 10. Wei, M. C.; Lindsten, T.; Mootha, V. K.; Weiler, S.; Gross, A.; Ashiya, M.; Thompson, C. B.; Korsmeyer, S. J. tBID, a membrane-targeted death ligand, oligomerizes BAK to release cytochrome c. Genes Dev. 2000, 14, 2060-2071. 11. Walensky, L. D.; Pitter, K.; Morash, J.; Oh, K. J.; Barbuto, S.; Fisher, J.; Smith, E.; Verdine, G. L.; Korsmeyer, S. J. A stapled BID BH3 helix directly binds and activates BAX. Mol. Cell 2006, 24, 199-210. 12. Desagher, S.; Osen-Sand, A.; Nichols, A.; Eskes, R.; Montessuit, S.; Lauper, S.; Maundrell, K.; Antonsson, B.; Martinou, J. C. Bid-induced conformational change of Bax is responsible for mitochondrial cytochrome c release during apoptosis. J. Cell Biol. 1999, 144, 891-901. 13. Edlich, F.; Banerjee, S.; Suzuki, M.; Cleland, M. M.; Arnoult, D.; Wang, C.; Neutzner, A.; Tjandra, N.; Youle, R. J. Bcl-x(L) retrotranslocates Bax from the mitochondria into the cytosol. Cell 2011, 145, 104116. 14. Eskes, R.; Desagher, S.; Antonsson, B.; Martinou, J. C. Bid induces the oligomerization and insertion of Bax into the outer mitochondrial membrane. Mol. Cell. Biol. 2000, 20, 929-935. 15. Goping, I. S.; Gross, A.; Lavoie, J. N.; Nguyen, M.; Jemmerson, R.; Roth, K.; Korsmeyer, S. J.; Shore, G. C. Regulated targeting of BAX to mitochondria. J. Cell Biol. 1998, 143, 207-215. 16. Wolter, K. G.; Hsu, Y. T.; Smith, C. L.; Nechushtan, A.; Xi, X. G.; Youle, R. J. Movement of Bax from the cytosol to mitochondria during apoptosis. J. Cell Biol. 1997, 139, 1281-1292. 17. Suzuki, M.; Youle, R. J.; Tjandra, N. Structure of Bax: coregulation of dimer formation and intracellular localization. Cell 2000, 103, 645-654. 18. Annis, M. G.; Soucie, E. L.; Dlugosz, P. J.; Cruz-Aguado, J. A.; Penn, L. Z.; Leber, B.; Andrews, D. W. Bax forms multispanning monomers that oligomerize to permeabilize membranes during apoptosis. EMBO J. 2005, 24, 2096-2103. 19. Gavathiotis, E.; Reyna, D. E.; Davis, M. L.; Bird, G. H.; Walensky, L. D. BH3-triggered structural reorganization drives the activation of proapoptotic BAX. Mol. Cell 2010, 40, 481-492. 20. Gross, A.; Jockel, J.; Wei, M. C.; Korsmeyer, S. J. Enforced dimerization of BAX results in its translocation, mitochondrial dysfunction and apoptosis. EMBO J. 1998, 17, 3878-3885. 36 ACS Paragon Plus Environment

Page 36 of 40

Page 37 of 40 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

Journal of Chemical Information and Modeling

21. Renault, T. T.; Dejean, L. M.; Manon, S. A brewing understanding of the regulation of Bax function by Bcl-xL and Bcl-2. Mech. Ageing Dev. 2017, 161, 201-210. 22. Czabotar, P. E.; Lessene, G.; Strasser, A.; Adams, J. M. Control of apoptosis by the BCL-2 protein family: implications for physiology and therapy. Nat. Rev. Mol. Cell Biol. 2014, 15, 49-63. 23. Scorrano, L.; Korsmeyer, S. J. Mechanisms of cytochrome c release by proapoptotic BCL-2 family members. Biochem. Biophys. Res. Commun. 2003, 304, 437-444. 24. Dai, H.; Smith, A.; Meng, X. W.; Schneider, P. A.; Pang, Y. P.; Kaufmann, S. H. Transient binding of an activator BH3 domain to the Bak BH3-binding groove initiates Bak oligomerization. J. Cell Biol. 2011, 194, 39-48. 25. Gavathiotis, E.; Suzuki, M.; Davis, M. L.; Pitter, K.; Bird, G. H.; Katz, S. G.; Tu, H. C.; Kim, H.; Cheng, E. H.; Tjandra, N.; Walensky, L. D. BAX activation is initiated at a novel interaction site. Nature 2008, 455, 1076-1081. 26. Czabotar, P. E.; Westphal, D.; Dewson, G.; Ma, S.; Hockings, C.; Fairlie, W. D.; Lee, E. F.; Yao, S.; Robin, A. Y.; Smith, B. J.; Huang, D. C.; Kluck, R. M.; Adams, J. M.; Colman, P. M. Bax crystal structures reveal how BH3 domains activate Bax and nucleate its oligomerization to induce apoptosis. Cell 2013, 152, 519-531. 27. Bleicken, S.; Classen, M.; Padmavathi, P. V.; Ishikawa, T.; Zeth, K.; Steinhoff, H. J.; Bordignon, E. Molecular details of Bax activation, oligomerization, and membrane insertion. J. Biol. Chem. 2010, 285, 6636-6647. 28. Dewson, G.; Kratina, T.; Czabotar, P.; Day, C. L.; Adams, J. M.; Kluck, R. M. Bak activation for apoptosis involves oligomerization of dimers via their alpha6 helices. Mol. Cell 2009, 36, 696-703. 29. Dewson, G.; Kratina, T.; Sim, H. W.; Puthalakath, H.; Adams, J. M.; Colman, P. M.; Kluck, R. M. To trigger apoptosis, Bak exposes its BH3 domain and homodimerizes via BH3:groove interactions. Mol. Cell 2008, 30, 369-380. 30. Dewson, G.; Ma, S.; Frederick, P.; Hockings, C.; Tan, I.; Kratina, T.; Kluck, R. M. Bax dimerizes via a symmetric BH3:groove interface during apoptosis. Cell Death Differ. 2012, 19, 661-670. 31. Oh, K. J.; Singh, P.; Lee, K.; Foss, K.; Lee, S.; Park, M.; Aluvila, S.; Kim, R. S.; Symersky, J.; Walters, D. E. Conformational changes in BAK, a pore-forming proapoptotic Bcl-2 family member, upon membrane insertion and direct evidence for the existence of BH3-BH3 contact interface in BAK homooligomers. J. Biol. Chem. 2010, 285, 28924-28937. 32. Gumbart, J.; Chipot, C.; Schulten, K. Free energy of nascent-chain folding in the translocon. J. Am. Chem. Soc. 2011, 133, 7602-7607. 33. Gumbart, J. C.; Chipot, C. Decrypting protein insertion through the translocon with free-energy calculations. Biochim. Biophys. Acta 2016, 1858, 1663-1671. 34. Henin, J.; Pohorille, A.; Chipot, C. Insights into the recognition and association of transmembrane alpha-helices. The free energy of alpha-helix dimerization in glycophorin A. J. Am. Chem. Soc. 2005, 127, 8478-8484. 35. Maity, A.; Sinha, S.; Ganguly, D.; Ghosh Dastidar, S. C-terminal tail insertion of Bcl-xL in membrane occurs via partial unfolding and refolding cycle associating microsolvation. Phys. Chem. Chem. Phys. 2016, 18, 24095-24105. 36. Singh, R. K.; Chamachi, N. G.; Chakrabarty, S.; Mukherjee, A. Mechanism of Unfolding of Human Prion Protein. J. Phys. Chem. B 2017, 121, 550-564. 37. Mancinelli, F.; Caraglia, M.; Budillon, A.; Abbruzzese, A.; Bismuto, E. Molecular dynamics simulation and automated docking of the pro-apoptotic Bax protein and its complex with a peptide designed from the Bax-binding domain of anti-apoptotic Ku70. J. Cell. Biochem. 2006, 99, 305-318. 38. Rosas-Trigueros, J. L.; Correa-Basurto, J.; Benitez-Cardoza, C. G.; Zamorano-Carrillo, A. Insights into the structural stability of Bax from molecular dynamics simulations at high temperatures. Protein Sci. 2011, 20, 2035-2046. 37 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

39. Jiang, Z.; Zhang, H.; Bockmann, R. A. Allostery in BAX protein activation. J. Biomol. Struct. Dyn. 2016, 34, 2469-2480. 40. Hsu, Y. T.; Youle, R. J. Nonionic detergents induce dimerization among members of the Bcl-2 family. J. Biol. Chem. 1997, 272, 13829-13834. 41. DeLano, W. L. The Pymol molecular graphics system, San Carlos, CA USA: DeLano Scientific: . 2002. 42. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. 43. Brooks, B. R.; Brooks, C. L., 3rd; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: the biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545-1614. 44. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781-1802. 45. Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E.; Mittal, J.; Feig, M.; Mackerell, A. D., Jr. Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone phi, psi and side-chain chi(1) and chi(2) dihedral angles. J. Chem. Theory Comput. 2012, 8, 3257-3273. 46. Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - an N.Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. 47. Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. Constant-Pressure Molecular-Dynamics Simulation - the Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613-4621. 48. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical-Integration of Cartesian Equations of Motion of a System with Constraints - Molecular-Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327341. 49. Fiorin, G.; Klein, M. L.; Henin, J. Using collective variables to drive molecular dynamics simulations. Mol. Phys. 2013, 111, 3345-3362. 50. Comer, J.; Gumbart, J. C.; Henin, J.; Lelievre, T.; Pohorille, A.; Chipot, C. The adaptive biasing force method: everything you always wanted to know but were afraid to ask. J. Phys. Chem. B 2015, 119, 1129-1151. 51. Maragakis, P.; Spichty, M.; Karplus, M. Optimal estimates of free energies from multistate nonequilibrium work data. Phys. Rev. Lett. 2006, 96, 100602. 52. Comer, J.; Phillips, J. C.; Schulten, K.; Chipot, C. Multiple-Replica Strategies for Free-Energy Calculations in NAMD: Multiple-Walker Adaptive Biasing Force and Walker Selection Rules. J. Chem. Theory Comput. 2014, 10, 5276-5285. 53. Roux, B. Statistical mechanical equilibrium theory of selective ion channels. Biophys. J. 1999, 77, 139-153. 54. Bakan, A.; Meireles, L. M.; Bahar, I. ProDy: protein dynamics inferred from theory and experiments. Bioinformatics 2011, 27, 1575-1577. 55. Tournier, A. L.; Smith, J. C. Principal components of the protein dynamical transition. Phys. Rev. Lett. 2003, 91, 208106. 56. Lou, H.; Cukier, R. I. Molecular dynamics of apo-adenylate kinase: a principal component analysis. J. Phys. Chem. B 2006, 110, 12796-12808. 57. Salmas, R. E.; Yurtsever, M.; Durdagi, S. Investigation of Inhibition Mechanism of Chemokine Receptor CCR5 by Micro-second Molecular Dynamics Simulations. Sci. Rep. 2015, 5, 13180. 58. Majumdar, S.; Ghosh Dastidar, S. Ligand Binding Swaps between Soft Internal Modes of alpha,beta-Tubulin and Alters Its Accessible Conformational Space. J. Phys. Chem. B. 2017, 121, 118-128.

38 ACS Paragon Plus Environment

Page 38 of 40

Page 39 of 40 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

Journal of Chemical Information and Modeling

59. Priya, P.; Maity, A.; Ghosh Dastidar, S. The long unstructured region of Bcl-xl modulates its structural dynamics. Proteins: Struct. Funct. Bioinform. 2017, 85, 1567-1579. 60. Newman, M. E. J.; Girvan, M. Finding and evaluating community structure in networks. Phys. Rev. E 2004, 69, 026113. 61. Sethi, A.; Eargle, J.; Black, A. A.; Luthey-Schulten, Z. Dynamical networks in tRNA:protein complexes. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 6620-6625. 62. Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graph. 1996, 14, 33-38. 63. Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562-12566. 64. Ovchinnikov, V.; Karplus, M.; Vanden-Eijnden, E. Free energy of conformational transition paths in biomolecules: the string method and its application to myosin VI. J. Chem. Phys. 2011, 134, 085103. 65. Pan, A. C.; Roux, B. Building Markov state models along pathways to determine free energies and rates of transitions. J. Chem. Phys. 2008, 129, 064107. 66. Wang, K.; Gross, A.; Waksman, G.; Korsmeyer, S. J. Mutagenesis of the BH3 domain of BAX identifies residues critical for dimerization and killing. Mol. Cell. Biol. 1998, 18, 6083-6089. 67. Ensing, B.; Laio, A.; Parrinello, M.; Klein, M. L. A recipe for the computation of the free energy barrier and the lowest free energy path of concerted reactions. J. Phys. Chem. B 2005, 109, 6676-6687. 68. Latzer, J.; Shen, T.; Wolynes, P. G. Conformational switching upon phosphorylation: a predictive framework based on energy landscape principles. Biochemistry 2008, 47, 2110-2122. 69. Johnson, Q. R.; Lindsay, R. J.; Nellas, R. B.; Fernandez, E. J.; Shen, T. Mapping allostery through computational glycine scanning and correlation analysis of residue-residue contacts. Biochemistry 2015, 54, 1534-1541. 70. Doshi, U.; Holliday, M. J.; Eisenmesser, E. Z.; Hamelberg, D. Dynamical network of residueresidue contacts reveals coupled allosteric effects in recognition, catalysis, and mutation. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 4735-4740. 71. Brinda, K. V.; Vishveshwara, S. A network representation of protein structures: implications for protein stability. Biophys. J. 2005, 89, 4159-4170. 72. Zhou, H.; Hou, Q.; Hansen, J. L.; Hsu, Y. T. Complete activation of Bax by a single site mutation. Oncogene 2007, 26, 7092-7102. 73. Gallivan, J. P.; Dougherty, D. A. Cation-pi interactions in structural biology. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 9459-9464. 74. Durrant, J. D.; de Oliveira, C. A.; McCammon, J. A. POVME: an algorithm for measuring bindingpocket volumes. J. Mol. Graph. Model. 2011, 29, 773-776.

39 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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

For Table of Contents Use Only

BIM

Binding

Remotely

BAX

Activation:

Regulates

Insights from the Free Energy Landscapes Souvik Sinha, Atanu Maity, Shubhra Ghosh Dastidar*

40 ACS Paragon Plus Environment

Page 40 of 40