PUMA, which belongs to the BH3-only protein family, is an intrinsically disordered protein (IDP). It binds to its cellular partner Mcl-1 through its BH3 motif, which folds upon binding into an α helix. We have applied a structure-based coarse-grained model, with an explicit Debye—Hückel charge model, to probe the importance of electrostatic interactions both in the early and the later stages of this model coupled folding and binding process. This model was carefully calibrated with the experimental data on helical content and affinity, and shown to be consistent with previously published experimental data on binding rate changes with respect to ionic strength. We find that intramolecular electrostatic interactions influence the unbound states of PUMA only marginally. Our results further suggest that intermolecular electrostatic interactions, and in particular non-native electrostatic interactions, are involved in formation of the initial encounter complex. We are able to reveal the binding mechanism in more detail than is possible using experimental data alone however, and in particular we uncover the role of non-native electrostatic interactions. We highlight the potential importance of such electrostatic interactions for describing the binding reactions of IDPs. Such approaches could be used to provide predictions for the results of mutational studies.
The discovery of intrinsically disordered proteins (IDPs) [1–3], requires a reconsideration of the principles of protein-protein interactions, which have largely been examined using well-defined structured proteins. IDPs are involved in many critical physiological processes, in particular within protein interaction networks such as transcriptional and translational regulation, cellular signal transduction, protein phosphorylation, and molecular assembly [4]. Whilst disordered/unstructured at physiological conditions, IDPs that play such roles often undergo conformational changes to upon binding to their biomacromolecular partners [4,5] in a process known as “coupled folding and binding” [5]. In the cell correct folding and binding of many proteins is essential for biological function. Since the 1890s the dominant descriptions of binding mechanisms are based on those described by Fischer: “lock and key” [6], “induced fit” [7] and “conformational selection” [8] (flexible binding). The lock and key mechanism, which involves rigid binding, is excluded for reactions that couple folding and binding. Recent kinetic studies have started to investigate the applicability of the “induced fit” and “conformational selection” mechanisms, which differ in whether the ligand is folded before or after binding. Electrostatic interactions are known to play important roles in binding processes [9–11]. However so far only a limited number of experiments have investigated these interactions within coupled folding and binding reactions [12–14]. One of the most extensively experimentally studied coupled folding and binding processes is that between PUMA and Mcl-1. It is already known that electrostatic interactions are involved in association of PUMA with Mcl-1 as association rates are affected by salt concentrations [12]. Here we study the process of coupled folding and binding of PUMA with Mcl-1 using coarse-grained simulations that are specifically designed to include contributions from electrostatic interactions to uncover the nature of their involvement. Mcl-1 is a stable folded pro-survival Bcl-2 protein that plays a critical role in development and tissue homeostasis. It is expressed in a range of tissue types, and is induced by a variety of stimuli to block apoptosis [15]. Gene knockout studies have demonstrated that Mcl-1 is required for embryonic and immune cell development [16,17], while its over-expression is implicated in cancer and resistance to cancer treatments [15]. Pro-apoptotic BH3-only proteins such as PUMA, many of which are intrinsically disordered in their unbound state [18], form α-helices upon binding to Mcl-1. Molecular dynamic simulations can gain many essential characteristics of IDPs that can not be obtained in experiments. However, the results of IDPs simulations are strongly influenced by the force field used [19–22]. Structure-based modeling can be a powerful tool for understanding coupled binding and folding of IDPs [23–27]. Classical structure-based models (SBMs) are established based on the native protein structure according to the conceptual framework of minimally frustrated energy landscapes. Therefore, the accuracy of IDP simulations is strongly dependent on parameters of SBM in the force field [27,28]. Simple SBMs, with only a single energy basin, may not be appropriate for the systems of IDP complexes or systems with many conformational states of similar energy. As a result, we applied an advanced SBM based on the stable IDP complex structure, and treated IDP and folded protein differently. For the accuracy of the simulations, we carefully recalibrated the strengths of intramolecular and intermolecular interactions based on the experimental structural characteristics of IDP in solution and the binding affinity of the IDP complex, which have been successfully applied in many IDP complex systems [25,27,29,30]. Besides, multi-basin models are also perfect choices to study the binding induced conformational changes of IDPs [31–34]. Here in this study, with a structure-based coarse-grained model [35–37] of the Mcl-1 and PUMA complex and a Debye—Hückel model to describe the electrostatic interactions at a moderate range of ionic strength, we investigate the contribution of ionic strength to the compactness of PUMA, and how this affects its rate of association with Mcl-1. Our results provide a residue level of detail for the mechanism of this modeled coupled folding and binding process, including the role of non-native electrostatic contacts in the formation of the encounter complex.
We aimed to describe the process of PUMA folding upon binding to its partner Mcl-1 in residue level detail using a simulation strategy that included consideration of electrostatic forces. To obtain the free energy landscape of the binding process between PUMA and Mcl-1, exchanges between unbound and bound states are needed in the thermodynamical simulations. Bias was added on the native contact pairs, aiming to enhance exchange rates between the unbound and bound states. After simulations, bias was removed to obtain the true free energy distribution. Firstly, the free energy of Mcl-1 and PUMA complex was projected on one reaction coordinate of intermolecular contact (Qinter) to provide an overview of the binding/unbinding process. A series of simulations were performed to calibrate the strength of intermolecular contacts (β). As shown in Fig 1B, when β = 0. 9, the binding free energy of Mcl-1·PUMA complex is converged to about 7. 93 kJ/mol, which matches the experimentally determined value of Kd [12,38]. Using β = 0. 9 we then investigated the nature of the coupled folding and binding process through two replicas of well-tempered metadynamics [39]. Fig 1A shows the free energy of the Mcl-1·PUMA complex as a function of the proportion of formed native intermolecular contacts (Qinter). There is a relatively small energy barrier of about 4 kJ/mol from the unbound state to the transition state, which is located at about 0. 1 of Qinter. In contrast the energy of the system declines very sharply towards the bound state, which is located at about 0. 72 Qinter. The bound state of Mcl-1·PUMA complex is so stable that the exchange between bound and unbound states is very slow (there are about only round 10 transitions from bound to unbound in each 128 ns simulation). Snapshots of PUMA binding to Mcl-1 with various Qinter are shown in Fig 1C to illustrate the process of binding although there are many potential structures for each value and simulations do not always involve the same contacts being formed first (see later). The free energy landscape for coupled folding and binding was further examined using additional reaction coordinates to examine the nature of the process. Taken together intermolecular and intramolecular contact formation (Qinter and Qintra) demonstrates a typical two-state system without populated intermediates (Fig 1) as has previously been suggested by experimental studies [12]. In its unbound state, the IDP PUMA is partly disordered with Qintra about 0. 4. Reflecting the coupling between folding and binding, within complex with Mcl-1 Qintra of PUMA increases to about 0. 7 (Qinter about 0. 72). Notably the Qintra of PUMA does not change considerably between the unbound and transition states, demonstrating that the majority of PUMA folding takes places after binding i. e. via an induced fit mechanism. We previously postulated this mechanism after a Φ-value analysis for PUMA binding to Mcl-1 showed little structure formation was present in the transition state [30]. However as noted in the article experimental kinetic measurements are not able to determine mechanism unambiguously as we are able to here. With the aim of maximizing comparison of our results with those obtained previously through experiment we also considered an further reaction coordinate; the helical content of PUMA. In the bound state the IDP PUMA is folded as a part of regular alpha helix so Qintra is related to the change of secondary structure of PUMA. Within the unbound and bound states the minima of the alpha helix content are about 10-20% and 40-60%, respectively. These results are consistent with the CD results of experiment [30]. In addition we are able to estimate the alpha helical content within the (unpopulated) transition state as approximately 18%. IDPs are typically enriched in charged residues, and the PUMA BH3 domain is no exception. Within our 34 amino acid peptide there are 10 negatively charged residues and 5 positively charged residues (PUMA sequence is included in supplementary materials S1 Fig). Such an abundance of charged residues suggests that electrostatic interactions may play an important role in determining the structure of PUMA in its unbound ensemble. We therefore examined the effect of increasing the ionic strength upon the structure of unbound PUMA. High salt concentrations act to shield charges from each other and hence reduce intramolecular electrostatic interactions. We calculated the radius of gyration (Rg) as a crude descriptor of the random coil ensemble of IDP PUMA. Fig 2 shows that Rg of PUMA does decrease with ionic strength, although the decrease is only small (from about 1. 34 nm to 1. 32 nm). Thus electrostatic interactions cause PUMA to be marginally more expanded than it would otherwise be. To probe this finding further we examined inter-residue distances within unbound PUMA. Distance maps (Fig 2B) reveal a lack of overall structural features. Instead inter-residue distance appears to be determined simply by separation distance in the amino acid sequence. This is highly consistent with its description as an IDP. However a more detailed examination demonstrates that electrostatic interactions modulate the inter-residue distances. Most notably, the distance maps for simulations with and without electrostatic interactions show that the distance between the two termini of PUMA increases slightly when charges are switched on. When electrostatic interactions are marginally shielded (IS = 10 mM) the distance between the two ends is about 2 Å larger when compared with no electrostatic interactions (Fig 2B). Other features are also apparent in the difference distance map, particularly for the arginine-containing region 12-18, which moves (up to) 1 Å closer to the negatively charged N-terminal region, and the region 18-30 when charges are introduced to the simulation. We further probed the simulations for evidence of more local changes by examining residue-specific electrostatic interactions within unbound PUMA. The distance between head and tail of PUMA is very large, suggesting a long extended coil with no hairpins/turns (as shown in Fig 2B). The change of structure is largely related to the local interactions. We considered contacts made between oppositely and similarly charged residues i. e. both attractive and repulsive intra-peptide interactions. For structures, we considered intra-PUMA contacts as fully formed (contact number 1) if the distance between the “atoms” was under 4. 5 Å. We considered intra-contacts as partially formed (contact number 0. 5) if the distance was larger than 4. 5 Å but under 6. 0 Å. The distribution of contact numbers for the two interaction types shows that considerably more repulsive contacts (8 on average) are formed than attractive contacts (4 on average) (Fig 3A). However most of the contacts are between directly connected charged residues (black, Fig 3C) that must always be within 6. 0 Å of each other. Excluding these there is an average of only one attractive and one repulsive contact being formed at any point (Fig 3B). Without electrostatic interactions, the probability of former opposite-charged contacts decreases, but the probability of former same-charged contacts increases. By contrast, the probability of intra-PUMA contacts when charges are switched off is shown in Fig 3A and 3B. A contact map showing the probability of each individual interaction being present demonstrates that most of the charge contacts are between residues relatively close together within the PUMA sequence (Fig 3C and S2 Fig). To reveal the role of electrostatic interactions on the binding behaviour of Mcl-1·PUMA complex, we performed kinetic simulations without explicit electrostatic interactions, and in the presence of electrostatic interactions at 6 different ionic strengths (from 10 to 250 mM). For each ionic strength 200 simulations with different initial velocities and structures were performed to ensure efficient sampling. The nature of our simulation approach meant that each simulation contained a single binding event. In the absence of electrostatic interactions this binding event took place on average 1330 ± 90 ps into each simulation, which is equivalent to a binding rate of (1. 3 ± 0. 1) × 109 s−1. As shown in Fig 4 this rate is increased when electrostatic interactions are included, and in an ionic strength dependent manner. The mean time increases quickly between 10 and 50 mM ionic strength and plateaus at higher ionic strengths, approaching the time observed in the absence of electrostatic interactions. At 10 mM ionic strength the binding takes place around 2. 0 times faster with a rate of (2. 6 ± 0. 1) × 109 s−1. To probe the basis of this phenomenon we then considered the binding process as two parts: an encounter (or capture) step and a further evolution step [41]. Broadly speaking binding rates could be modulated through changes in three separate key parameters: rate of formation of encounter complex (capture), rate of evolution from encounter complex to bound complex, and probability of progression from encounter complex to bound complex. We examined the dependence of each of these separate stages on the ionic strength. PUMA formed multiple encounter complexes before successful binding occurred. The rate of capture events without electrostatic interactions is (2. 2 ± 0. 1) × 1010 s−1. Similarly to the binding rate, this is increased in the presence of electrostatic interactions in an ionic strength dependent fashion, with a 2. 1-fold rate enhancement at 10 mM ionic strength, (4. 6 ± 0. 2) × 1010 s−1 (Fig 4). Thus the binding rate and capture rate have remarkably similar behaviour with ionic strength, suggesting that electrostatic forces enhance binding rates essentially entirely through altering collision rates. Consistent with this proportion of successful capture events appears independent of ionic strength (Fig 4D). The evolution rate is also relatively insensitive to electrostatic interactions, ranging only from (4. 3 ± 0. 2) × 1010 s−1 without electrostatic interactions to (3. 6 ± 0. 14) × 1010 s−1 at 10 mM ionic strength (Fig 4C). For PUMA systems at low ionic strength (ionic strength = 10 mM), we have shown electrostatic interactions can accelerate the binding process. However there are very few native electrostatic intermolecular contacts within the encounter complex in our simulations (average number of contacts is only 0. 1). However, the average of native charged contacts within 6 Å is 0 (Fig 5B), which suggests that the distances of the native charged contacts are all longer than 6. 0 Å. The contact map of native opposite-charged contacts (Fig 5D) indicates that there is only one native opposite-charged contact with high probability, which is between residue Lys66 of Mcl-1 and Glu9 of PUMA (Fig 6A). In contrast there is a much wider distribution for the number of non-native opposite-charged contacts, with a much higher average contact number (5. 1) than that for native contacts. This contact number is also somewhat higher than that for non-native same-charged contacts (3. 5) (Fig 5A). The non-native charged interactions are almost all long-range (6. 0—10. 0 Å) (Fig 5A and 5B), and whilst it appears that many different interactions can be present in the encounter complex, some interactions are more likely than others. In particular interactions between the negatively charged residues in the N-terminal region of PUMA and Arg65/Lys66 of Mcl-1 are all present in over 20% of encounter complexes, and interactions between Arg15/Arg16 of PUMA and Glu57 of Mcl-1, and between Glu3/Glu4 of PUMA and Arg80 of Mcl-1 are present in over 10% of encounter complexes (see Figs 5D and 6A). Interestingly the N-terminal region of PUMA is also the region which we were previously able to show had formed some structure within the transition state for binding [30]. Therefore, PUMA firstly gets close to Mcl-1 with long-ranged electrostatic formed by its N-terminal residues (negatively charged) with positively charged residues (Arg65, Lys66, and Arg80) of Mcl-1. There are more native non-electrostatic contacts (2. 1) than native electrostatic contacts (0. 1) in the encounter complex. Likewise, all the native non-electrostatic contacts are long-ranged (distance longer than 6. 0 Å, Fig 5C). We note that there are also many non-native non-electrostatic interactions within the encounter complex (average 7. 5, Fig 5C). However, none of the non-native non-electrostatic interactions has a probability > 0. 05 (S3 Fig). To dissect more fully the role of electrostatic forces from non-electrostatic ones in formation of the encounter complex we also performed simulations in the absence of charges for comparison. As previously very few native contacts were formed between previously charged residues. However the average number of non-native electrostatic contacts, both “opposite-charged” and “same-charged”, decreased to 2. 1 and 1. 9 respectively. The probability distribution of contact numbers has also changed so that the most likely number of contacts is 0 (Fig 7A and 7B), rather than 7 and 4 as obtained in simulations including electrostatics (10 mM). Most of the contacts formed are long-range. Examining the changes at the residue specific level indicates that this is mainly a result of a large decrease in “opposite-charged” contacts between Glu2/Glu3/Glu4/Glu5 of of PUMA and Mcl-1 (Fig 7D and S4A Fig), which will also affect the non-electrostatic interactions at this region (S4C Fig).
We used coarse-grained molecular dynamics simulations to describe the mechanism of interactions between the natural intrinsically disordered PUMA BH3 peptide and its Bcl-2 partner protein Mcl-1. Our findings portray an induced fit binding process, with PUMA being relatively disordered at the transition state. This is in good agreement with previous experimental kinetic analyses that have suggested a relative unstructured transition state on the basis that the dissociation rate constant is much more sensitive to denaturant concentration than the association rate constant is [12]. Our simulations suggest the helical content of the transition state is approximately 18%, much closer to that of minima of the unbound state (10-20%) than the bound state (40-60%). This is highly consistent with, but extends, our knowledge of helical contents on the folding landscape from experimental studies that are fundamentally limited to populated states [30]. This compares favourably with Leffler α = 0. 10 ± 0. 01 determined from a recent Φ-value analysis, which aims to estimate the proportion of global structure formation at the transition state [30]. This work also showed the transition state to contain almost no structure at the C-terminus and only partial structure at the N-terminus [30]. This is highly consistent with the view formed of the encounter complex, formed before the transition state, from our simulations. These show the most probable contacts are formed between the N-terminus of PUMA with Mcl-1 residues that are close to its final position in the bound state. Inter-residue distance maps of PUMA support its overall structural characterisation as an IDP since the largest determinant of the distance between any two residues is their sequence separation. Chebaro et al. [42] obtained a description of the energy landscape of PUMA and showed that the helix was unstable in the unbound state. The high proportion of charged residues within the PUMA sequence prompted us to examine whether structural changes might result from altering the ionic strength of the surrounding medium. Within our sequence of PUMA, there are 66 potential opposite-charged contact pairs and 70 potential same-charged contact pairs (including termini). The fraction of positive charged residues (f+) is 0. 147 and the fraction of negative charged residues (f−) is 0. 294. Then, the fraction of charged residues (FCR) and the net charge per residue (NCPR) /mean net charge are 0. 441 and 0. 147 (FCR = (f+ + f−); NCPR = |f+ − f−|), respectively. The mean hydrophobicity of PUMA is computed via the online tools ProtScale of ExPASy [43] (http: //web. expasy. org/protscale/) with a window size of 5 residues and rescaled to fit between 0 and 1, which is 0. 340. In addition, the intra-chain distance (Rij) has a relationship with chain separation (|i − j|), which can fit to the line of y = 0. 331 x0. 683. As shown in the Uversky diagram [44] and Das and Pappu diagram [45] of supplementary materials S5 Fig, PUMA is classified as the strong polyampholytes of IDPs [45–47]. Therefore, the characteristic of PUMA sequence is consistent with its behavior of Rg at different ionic strengths. However both measures of large-scale structural tendencies we employed, Rg and inter-residue distance maps, demonstrated only relatively small changes due to electrostatic interactions. Notably after 150 mM ionic strength there is no further decrease in radius, so the physiological situations may be well represented by simulations that do not include electrostatic contributions explicitly. We are able to examine the subtle changes in PUMA structure in more detail. Unbounded PUMA is shown to be a extended coil. However the most significant contribution to the minor increase in PUMA radius when electrostatics are included in simulations, appears to be interactions between the two ends of the PUMA peptide. Both termini are negatively charged overall due to the presence of multiple glutamic acid residues (Glu2, Glu3, Glu4, Glu5, Glu31 and Glu32), and the difference inter-residue distance map indicates these move around 2 Å apart when charges are introduced and only slightly screened (IS = 10 mM). Indeed, the C-terminus (residue 30-34) moves further apart from the rest of the peptide in general. We have shown that intra-chain electrostatic interactions within PUMA enlarge its radius. Thus at low ionic strengths the radius of interaction between Mcl-1 and PUMA may be increased. For induced fit processes such increased capture radii has been theoretically proposed to accelerate binding processes, a phenomenon termed the “fly-casting effect” [48,49]. This phenomenon is thought to alter association rates less than 2-fold when comparing folded and random coil versions of the same protein. Thus small modifications in the radius such as those observed here may be expected to produce only tiny changes in association rates. Nonetheless we note that particular structural changes, especially those that alter the accessibility of charged residues, could impact on the efficiency of any “steering” of a protein towards its partner. Electrostatic interactions lead to both attractive and repulsive forces. Our simulations indicate that in the case of the PUMA⋅Mcl-1 complex the binding process is accelerated by electrostatic interactions i. e. attractive forces dominate. We observed that association rates decreased at higher ionic strengths where charges are more efficiently screened from each other and electrostatic interactions weakened. This behaviour is similar to that observed in many kinetic experiments of protein-protein association where association rates are found to vary depending upon the salt concentration. In fact such rate enhancement has already been observed experimentally for PUMA association with Mcl-1, though the effect was slightly more pronounced [12]. Whilst in experiment kon decreases roughly 2-fold from 10 mM to 50 mM (about 9. 8 × 107 and 4. 4 × 107 M−1 s−1), in our simulations, the binding rate decreases only just over 1. 3-fold over the same range. This may be due to the peptide differing by a single charged residue from that used in the experiments, or missing fine structural details as a result of our coarse-grained molecular dynamics approach. Examining this behaviour through simulation allows us to separately consider the role of electrostatic interactions in capture and evolution parts. Our data suggest that electrostatic forces assist the initial formation of the encounter complex, and play a lesser role in subsequent progression to the bound state. Simulations from the group of Jianhan Chen for p27 with Cdk2/cyclin A [26], and for p53-TAD1 with TAZ2, HIF-1a with TAZ1, and NCBD with ACTR [50] all demonstrated that electrostatic forces act to reduce energy barriers between the various states, thus enhancing collision rates. In contrast to our findings they also suggested an important role of electrostatic contacts in progression from the collision complex to the fully bound state [26,50]. In fact we observe a 1. 2-fold decrease in evolution rate upon introducing partially screened charges to the simulation, suggesting electrostatic interactions are actually retarding progression from the encounter complex to the bound state. Nonetheless this change remains smaller than the effect upon binding and capture rates, and given that the evolution time remains appreciably lower than the timescale for binding this is unlikely to have a meaningful contribution to binding rates in our simulations. Since the effect of altering ionic strength is mainly on the long-range electrostatic interactions of capture process we investigated the nature of the charge interactions within the encounter complex. In our model of the Mcl-1·PUMA complex there are only 3 native opposite-charged contacts (between Lys66 of Mcl-1 and Glu9 of PUMA, Asp88 of Mcl-1 and Arg15 of PUMA, Arg95 of Mcl-1 and Asp19 of PUMA) and 3 native same-charged contacts (between Lys87 of Mcl-1 and Arg15 of PUMA, Asp88 of Mcl-1 and Asp19 of PUMA, Arg95 of Mcl-1 and Arg15 of PUMA). Rate constants for binding reactions between pairs of folded proteins have been described to depend strongly upon charges located in the binding interface [51]. As an intrinsically disordered protein that undergoes induced fit binding to Mcl-1, PUMA does not have a well-defined binding interface in its unbound state. Nonetheless if the native interfacial contacts are involved in electrostatic rate enhancement of the capture rate then they will be present within the encounter complex. In contrast we observed only one native opposite-charged interaction (between Lys66 of Mcl-1 and Lys9 of PUMA) with long-range distance (distance larger than 6 Å) in encounter complex, which indicates the entrance for PUMA. Most native interactions formed in encounter complex are non-electrostatic in nature. Long-range electrostatic interactions do play an important role in the capture event though, as evidenced by the ionic strength dependence of the capture rate. Instead non-native opposite-charged interactions contribute a lot on the capture process. On average over 5 such contacts are formed in each encounter complex at 10 mM ionic strength. There is a wide distribution of probable contact numbers, ranging from 1 to 7, which reflects the non-uniformity of encounter complex structures. The most commonly observed opposite-charged contacts are between Arg65/Lys66/Arg80 of Mcl-1 and Glu2/Glu3/Glu4/Gl5 of PUMA. Each charged residue has been found to form multiple different intermolecular contacts within the encounter complex, however mapping these contacts onto the final structure demonstrates that contacts are more likely where residues are close to each other in the final structure. Similar coarse-grained simulations of the Mcl-1 and PUMA complex have been performed earlier by Rogers et al. [30]. The N-terminus was demonstrated to be important for the early stages of PUMA binding, even without electrostatic interactions. Here we show further that the non-native electrostatic interactions between N-terminus of PUMA and Mcl-1 can act as a “steering force” for binding, especially at low ionic strengths. A similar finding was reported recently for a peptide of RA-GEF2 binding to its PDZ domain partner where various non-native dynamic salt bridges were found within the encounter complex [52]. Further evidence for the importance of these contacts is that excluding electrostatic forces from our simulations alters the distribution of PUMA over the Mcl-1 surface in the encounter complex, with less focussed “charged interactions”.
NMR coordinates of Mcl-1 complexed with PUMA (PDB ID 2ROC [53]) were prepared for constructing the model. In this structure, the PUMA is 27 a. a. length. According to the experiments of Rogers et al. [12], a full-length PUMA (34 a. a.) was build by using Chimera software [54] (S1 Fig). The initial coarse-grained Cα structure-based model (SBM) of Mcl-1 and PUMA complex was generated using SMOG on-line toolkit, which included one bead on the Cα atom of each residue of the complex [35,36,55,56]. This model contains Mcl-1 (162 residues) and PUMA (34 residues), 424 and 36 intramolecular contacts within each of them, as well as 77 intermolecular contacts between them. The native contact map was built by the Shadow Algorithm [55]. The potential energy function consists of both bonded and nonbonded terms. Additionally, we introduced the charge characterization into our SBM model to study the electrostatic interactions in this system. As a result, the potential energy form used in this study is given in the following equation: V = ∑ b o n d s ϵ r (r - r 0) 2 + ∑ a n g l e s ϵ θ (θ - θ 0) 2 + ∑ d i h e d r a l K ϕ (n) 1 - cos (n × (ϕ - ϕ 0) ) + ∑ c o n t a c t s ϵ i j 5 σ i j r i j 12 - 6 σ i j r i j 10 + ∑ n o n - c o n t a c t s ϵ N C σ N C r i j 12 + V D e b y e - H u ¨ c k e l (1) In Eq 1, ϵr = 100 ϵ, ϵθ = 20 ϵ, K ϕ (1) = ϵ and K ϕ (3) = 0. 5 ϵ. The interaction strength of Lennard-Jones-type potential is proportional to the statistical potential reported for the residue types of i and j by Miyazawa and Jernigan to generate the flavored model. [57] Therefore, the coefficient of nonbonded contacts, ϵij is set as follows: ϵ i j = γ ϵ i j M J ϵ ¯ M J - 1 + 1 (2) where ϵ i j M J is the original MJ potential, ϵ ¯ M J is the mean value of the entire set of MJ weights in the complex system, and γ is set to 1. 0 corresponding to the “flavored model” [58]. As with similar strategies in previous studies [25,27,29], the native nonbonded potential can be separated into intramolecular and intermolecular terms, which should be rescaled according to the experimental data (see the part of Method in supplementary materials). The electrostatic interaction is calculated by the Debye—Hückel model, which can quantify the strength of charge-charge attractions and repulsions at various salt concentrations: V D e b y e - H u ¨ c k e l = Γ D H × K c o u l o m b B (κ) ∑ i, j q i q j exp - κ r i j ϵ r i j (3) In Eq 3, Kcoulomb = 4 πϵ0 = 138. 94 kJ·mol−1·nm·e−2 is the electric conversion factor; B (κ) is the salt-dependent coefficient; κ−1 is the Debye screening length, which is directly influenced by the solvent ionic strength (IS) /salt concentration Csalt (κ ≈ 3. 2 C s a l t); ϵ is dielectric constant, which is set to 80 during the simulations. ΓDH is the energy scaled coefficient which aims to make the total energy balanceable. In our model, Lys and Arg have a positive point charge (+e), Asp and Glu have a negative point charge (-e). All the charges are placed on the Cα atoms. Besides the systems in altering ionic strengths, there is also a system with no electrostatic interactions. Under physiological ionic strengths (Csalt ∼ 0. 15 M), κ is 1. 24 nm−1, so we set ΓDH = 0. 535 in our simulations, so that VDH for two opposite charged atoms located at a distance of 0. 5 nm matches the native contact energy. More details of Debye—Hückel model can be found in these papers [34,59–61]. Our definition of native contact cutoff is the same as that used in many previous structure-based model simulations, and is based on the distances between the residues in the experimental structures. In contrast, there is no universal standard for the definition of the cutoffs of non-native contacts. It is important that we analyze the native and non-native interactions in a consistent manner. For intermolecular non-native contacts, we checked all the native distances of the native contacts. We find that in the native complex the average distances between two Cα atoms in two chains that form a contact is about 8. 6 Å. In addition, the minimum of native distance between two Cα atoms in two chains is 4. 9 Å. As a result, when analysing contact number within the encounter complex we considered a contact between Mcl-1 and PUMA as fully formed (contact number 1) if the distance between the “atoms” was less than or equal to 6. 0 Å (about 1. 2 times of lowest distance of native inter-chain contacts); we considered contacts as partially formed (contact number 0. 5) if the distance was larger than 6. 0 Å but under 10. 0 Å (non-native contacts, about 1. 2 times of mean distance of native inter-chain contacts) [38]. This method was chosen to allow the same way of counting contact number for both native and non-native interactions. ln (bind rate), ln (capture rate), and ln (evolution rate) were obtained from FPTon, MPTcap, and FPTevo, respectively, ln k = 1 n ∑ i n ln 1 t a u i (4) here, k represents bind rate/capture rate/evolution rate, tau represents FPTon/MPTcap/FPTevo, and n is 200 runs at each ionic strength. Collision success probability is obtained from bind rate and capture rate, c o l l i s i o n s u c c e s s p r o b a b i l i t y = b i n d r a t e c a p t u r e r a t e I b i n d r a t e c a p t u r e r a t e 0 (5) where I is the ionic strength.
Intrinsically disordered proteins (IDPs) are typically enriched in charged residues, with a general “coupled folding and binding” mechanism. We uncover the importance of non-native long-range electrostatic interactions between PUMA and Mcl-1 in the early stage of binding, by using the coarse-grained molecular dynamic simulations with charged model, which reveals the reason of the obviously increased binding rate at low ionic strength. However, intramolecular electrostatic interactions within PUMA can not influence the folding extent of unbound states of PUMA significantly. The results of binding-folding energy landscapes support an “induced fit” mechanism of binding. The critical residues and interactions between PUMA and Mcl-1 suggested here will be helpful for further mutational and mechanism studies.
ionic binding method
chemical bonding
molecular dynamics
electrostatic bonding
electricity
electrostatics
molecular biology techniques
thermodynamics
intrinsically disordered proteins
reaction dynamics
physical chemistry
research and analysis methods
enzyme immobilization
proteins
chemistry
molecular biology
free energy
physics
biochemistry
biochemical simulations
computational chemistry
transition state
biology and life sciences
physical sciences
computational biology
Role of non-native electrostatic interactions in the coupled folding and binding of PUMA with Mcl-1