Comparing effects of attractive interactions in crowded systems: nonspecific, hydrophobic, and hydrogen bond interactions

Coarse-grained model for protein folding

To model protein folding, we use the coarse-grained “C_β-model” described in Bhattacherjee & Wallin (2012). It is a model with three different types of amino acids: polar (p), hydrophobic (h), and turn (t) amino acids. Geometrically, the protein chain is described using an atomistic backbone (C_α, C′, N, H, H_α1, and O) and simplified sidechains using an enlarged C_β atom. The t amino acids differ from p and h in that it does not contain a C_β atom, which is instead replaced by an H_α2 atom. Hence, t is strongly related to glycine and more flexible than p and h. Chain conformations are completely specified by the 2N Ramachandran angles ${\left\{{\varphi }_{\text{i}},{\psi }_{\text{i}}\right\}}_{\text{i}=1}^N$. Hence, bond lengths, bond angles, and dihedral angles (e.g., the peptide plane angle ω = 180°) are held fixed at standard values. The potential energy function of the model, E_p, includes four terms E_p = E_loc + E_exvol + E_hbond + E_hp, described in detail in Bhattacherjee & Wallin (2012). Briefly, the first term, E_exvol, provides atom-atom repulsions with a range σ_a determined by the sum of the atoms’ van der Waals radii, ${\sigma }_{\text{a}}={\sigma }_{\text{i}}^{\text{vdW}}+{\sigma }_{\text{j}}^{\text{vdW}}$. The second term, E_loc, represents electrostatic interactions between partial charges of adjacent peptide planes. This term is included to ensure sampled ϕ_i, ψ_j angles agree with statistics from real protein structures (Ramachandran plots).

The two last terms, E_hbond and E_hp, represent hydrogen bonding and effective hydrophobic interactions, respectively. These terms are essential for stabilizing and driving the formation of the native state. Hydrogen bonds are modeled using: (1)${\displaystyle E_{\text{hbond}}=k_{\text{hb}}\sum _{\text{ij}}{\gamma }_{\text{ij}}\left[5{\left(\frac{{\sigma }_{\text{hb}}}{r_{\text{ij}}}\right)}^{12}-6{\left(\frac{{\sigma }_{\text{hb}}}{r_{\text{ij}}}\right)}^{10}\right]\times \left(cos{\alpha }_{\text{ij}}cos{\beta }_{\text{ij}}\right)\frac{1}{2},}$

where the sum is over pairs of CO, NH groups, r_ij is the OH distance, σ_hb = 2.0 Å, and k_hbond = 3.1. The interaction strength is modified by a sequence-dependent scale factor γ_ij taken to be 1 for hh, hp, and pp pairs, and 0.75 for tt, th, and tp pairs. The reduced hydrogen bonding capacity of t amino acids mimics the secondary structure breaking ability of glycine. The directional dependence is implemented via the factor ${\left(cos{\alpha }_{\text{ij}}cos{\beta }_{\text{ij}}\right)}^{\frac{1}{2}}$, where α_ij and β_ij are the NHO and HOC′ angles, respectively. Additionally, if either α_ij < 90°or β_ij < 90°, the ij contribution is set to zero. The hydrophobic effect is modeled using pairwise additive interactions between hydrophobic amino acids. Specifically, (2)${\displaystyle E_{\text{hp}}=-k_{\text{hp}}\sum _{\text{ij}}g\left(r_{\text{ij}};{\sigma }_{\text{hp}}\right),}$ where the sum is over pairs of C_β atoms of h amino acids, excluding nearest and next-nearest neighbors along the chain, and (3)${\displaystyle g\left(r;r^0\right)=exp\left(-\frac{{\left(r-r^0\right)}^2}{2}\right)}$ is a function with a maximum at r = r⁰. We set σ_hp = 5 Å and the interaction strength k_hp = 0.805. In this model, tuning the relative strength of hydrogen bonding and hydrophobic forces is essential to obtain folding into stable and protein-like structures (Irbäck, Sjunnesson & Wallin, 2000; Irbäck, Sjunnesson & Wallin, 2001).

Crowders

Repulsive interactions are modeled using the pair potential (Mittal & Best, 2010) (4)${\displaystyle V\left(r\right)={\left(\frac{\sigma }{r-\rho +\sigma }\right)}^{12},}$ where r is a crowder-crowder or crowder-atom distance. Because V(ρ) = 1 and V → ∞ as r → ρ − σ (for r < ρ − σ, we set V = ∞) the two parameters ρ and σ determine the range and “sharpness” of the repulsion, respectively. We set ρ = 2R_c for crowder-crowder interactions and ρ = ρ_cp = R_c + σ_a for crowder-atom interactions, where R_c and σ_a are the radii of the crowder and atom, respectively. We set σ = 6 Å for crowder-crowder interactions and 3 Å + σ_a for crowder-atom interactions. A crowder-atom attractive interaction is modeled with the potential V(r) − ϵ_attg(r; ρ_cp), where the function g is given by Eq. (3). The form of this potential for different attraction strengths ϵ_att is shown in Fig. 1. For a system consisting of a protein with N_a atoms and N_cr crowder particles, the total potential energy then becomes E = E_p + E_cc + E_cp, where (5)${\displaystyle E_{\text{cc}}=\sum _{\text{i}}^{N_{\text{cr}}}\sum _{\text{j=i+1}}^{N_{\text{cr}}}V\left(r_{\text{ij}}\right)}$ and (6)${\displaystyle E_{\text{cp}}=\sum _{\text{i}}^{N_{\text{cr}}}\left[\sum _{\text{j}}^{N_{\text{a}}}V\left(r_{\text{ij}}\right)-{ϵ}_{\text{att}}\sum _{\text{j}}g\left(r_{\text{ij}};{\rho }_{\text{cp}}\right)\right]}$ are the crowder-crowder and crowder-protein energies, respectively. In Eq. (6), the second sum within square brackets controls which atoms are subject to crowder attraction. For crowders with nonspecific attraction to the protein, the sum is over all C_β atoms (p and h amino acids). For hydrophobic crowder, the sum is over the C_β atoms of h amino acids. For excluded volume crowders, ϵ_att = 0.

Simulated tempering Monte Carlo

To find the thermodynamic behavior of various protein-crowder systems, as determined by the amino acid sequence, number of crowders, and the energy function E(r), we use simulated tempering Monte Carlo (MC) (Marinari & Parisi, 1992; Lyubartsev et al., 1992; Irbäck & Potthast, 1995). In addition to a random walk in conformational space, as in basic Monte Carlo, simulated tempering also carries out a random walk in temperature while keeping the simulation at equilibrium. This is achieved by defining a set of temperatures, ${\left\{T_{\text{j}}\right\}}_{\text{j}=1}^{\text{M}}$, and simulating the joint probability distribution (7)${\displaystyle P\left(r,j\right)\propto e^{-{\beta }_{j}E\left(r\right)+g_j},}$ where β_j = 1/k_BT_j, k_B is Boltzmann’s constant, and j has been made a dynamic parameter. The g_j’s are M simulation parameters that control the marginal distribution P(j). Jumps between temperatures, j → j′, are accomplished as MC updates, with acceptance probability (8)${\displaystyle P_{\text{acc}}\left(r,j\to j^{\prime }\right)=min\left[1,e^{-E\left(r\right)\left({\beta }_{j^{\prime }}-{\beta }_j\right)+g_{j^{\prime }}-g_j}\right].}$ A common choice for the g_j parameters, which we follow here, is to select them such that P(j) is roughly flat, ensuring that sampling of conformations takes place at each temperature T_j.

Simulations and analysis details

Equilibrium behaviors of crowder-peptide systems are determined using simulated tempering Monte Carlo simulations (Marinari & Parisi, 1992). Runs are carried out with either 8 temperatures in the range k_BT = 0.48–0.68 or 10 temperatures in the range 0.40–0.70. For each system, 10 independent runs with 5  × 10⁹ elementary MC updates are performed. Initial configurations are obtained by picking random chain conformations and random spherical crowder positions, followed by a relaxation step to remove any hard-core steric clashes. Monte Carlo chain updates are divided equally between the protein chain and the crowder particles. Two different types of chain moves are used: biased Gaussian steps (BGS) (Favrin, Irbäck & Sjunnesson, 2001), which produce approximately local chain deformations, and pivot moves, which produce global changes. In the latter type, a single ϕ_i or ψ_i angle is set to a new random value such that the chain rotates around the NC_α or C_αC′ bonds. The frequencies of chain updates are set so that BGS is most frequent at low Ts and pivot is most frequent at high Ts. Temperature updates are attempted every 100 MC step. Spherical crowder simulations are carried out in a cubic box of side length L = 100 Å while for simulations of polypeptide crowders L = 65 Å. The number of spherical crowding agents are 14, 28, 42, and 56, corresponding to packing fractions ϕ_c = 0.10, 0.20, 0.30, and 0,40, respectively. The crowder radius is R_c = 12 Å. For simulation carried out with 8 temperatures in the range k_BT = 0.48–0.68, the multistate Bennett acceptance ratio (MBAR) technique was applied to determine thermodynamic averages in the range k_BT = 0.40–0.70.

Observables

We quantify native state stability using two quantities, ΔF and T_m. The free energy of folding is determined using (9)${\displaystyle \Delta F=F_{\text{N}}-F_{\text{U}}=-k_{\text{B}}Tln\frac{P_{\text{nat}}}{1-P_{\text{nat}}},}$ where F_U and F_N are the free energies of the unfolded and native states, respectively, and P_nat is the population of the native state. The native state population P_nat is determined as in Ref. Bazmi & Wallin (2022). Conformations are considered part of the native state if Q_nat ≥ Q_cut where Q_nat is the number of formed native contacts and Q_cut = 50. The folding midpoint temperature, T_m, is determined by fitting the temperature dependence of P_nat to a two-state model with two fit parameters.

Spherical crowders

Using MC simulations, we determined the thermodynamics behavior of our model protein, α₃₅, in isolation and in the presence of three different types of spherical crowders having the same radius (R_c = 12 Å) but different interactions with the protein chain: (1) excluded volume crowders, capable of only steric repulsions; (2) “nonspecific” crowders, which in addition to steric repulsions, form energetically favorable contacts with polar and nonpolar amino acids; and (3) hydrophobic crowders, which differ from the nonspecific crowders in that they interact favorably with nonpolar amino acids only. We varied the crowder volume fraction, ϕ_c, and the strength ϵ_att of the attractive protein-crowder interactions. In general, we find that the equilibrium stability of the α₃₅ native state, as quantified by the free energy of folding, ΔF =  − k_BTln[P_nat/(1 − P_nat)], where P_nat is the population of the native state, depends on both ϕ_c and ϵ_att. As an illustration of our results, Fig. 2 shows the temperature dependence of ΔF at ϵ_att = 1.5, for different ϕ_c and crowder types.

In the absence of crowders (solid black curve in Fig. 2), the helical hairpin fold of α₃₅ is highly stable at low T. For example, at the lowest studied temperature (k_BT = 0.40), which corresponds to $T\approx 0.75T_{\text{m}}^0$, where $T_{\text{m}}^0$ is the midpoint temperature of the α₃₅ folding transition for ϕ_c = 0, ΔF/k_BT =  − 7.1, which is within the range of stabilities typically found for single-domain proteins (Zeldovich, Chen & Shakhnovich, 2007). The midpoint temperature $T_{\text{m}}^0$, i.e., the temperature for which ΔF = 0, is also commonly used as a measure of protein stability. For α₃₅, $k_{\text{B}}{T}_{\text{m}}^0=0.535$, as obtained previously (Bazmi & Wallin, 2022) by fitting the of α₃₅ folding curve i.e., P_nat as function of T, to a two-state model with two free parameters. For the values of ϵ_att considered here, α₃₅ remains stable in the low T region even when the crowder packing fraction reaches ϕ_c = 0.40. However, there are clear variations between different crowder types as can be seen in Fig. 2. For example, for ϕ_c = 0.30 and k_BT = 0 .40, the excluded volume crowders give ΔF/k_BT =  − 8.3, which is a decrease relative to the ϕ_c = 0 case. For the nonspecific and hydrophobic crowders the corresponding ΔF/k_BT values are −6.1 and −3.3, respectively, which, by contrast, are increases relative to ϕ_c = 0.

At a temperature just below the folding midpoint, $T^{-}=0.95T_{\text{m}}^0$, we observe the following trends. Upon the addition of excluded volume crowders, ΔF decreases monotonically with ϕ_c, as shown in Figs. 3A and 3C, meaning a strict stabilization of the protein. The soft interactions of the nonspecific and hydrophobic crowders are expected to be destabilizing. However, because these two crowder types occupy the same space as the excluded volume crowders, they also provide a stabilizing effect. The net effect will be determined by a competition between the stabilizing steric repulsions and the destabilizing soft interactions. Indeed, at low attraction strength, ϵ_att = 1.0, the addition of either nonspecific or hydrophobic crowders leads to an overall stabilization of α₃₅ (see Fig. 3A and 3C). At ϵ_att = 1.5, the hydrophobic crowders overcome the stabilizing excluded volume effect leading to a net destabilization, while the nonspecific crowders are still net stabilizing (see Figs. 3B and 3D), although for the nonspecific crowders ΔF ≈ 0 at ϕ_c = 0.40. The picture obtained is similar if instead T_m is used to quantify native state stability. As shown in Fig. 3C, there is a monotonic increase in T_m for both the excluded volume crowders and for the hydrophobic and nonspecific crowders with ϵ_att = 1.0. For ϵ_att = 1.5, the nonspecific crowders still exhibit a monotonic increase in T_m with increases ϕ_c while the hydrophobic crowders instead exhibit a monotonic decrease (see Fig. 3D). Overall, our results show that excluded volume crowders provide a stabilizing effect on α₃₅, which can be counteracted by the presence of soft attractive interactions. The hydrophobic crowders provide a stronger destabilizing effect than the nonspecific crowders.

How can the stronger destabilizing effect of hydrophobic attractions be explained?

To address the question of why hydrophobic attractive interactions are more destabilizing than nonspecific interactions, we consider the interaction energy between crowders and protein, E_cp (see Methods). E_cp is a mix of repulsive and attractive energy contributions. We consider, in particular, $E_{\text{cp}}^{\text{U}}$ and $E_{\text{cp}}^{\text{N}}$, i.e., the crowder-protein energy determined separately for the U and N states. When the difference $\Delta E_{\text{cp}}=E_{\text{cp}}^{\text{U}}-E_{\text{cp}}^{\text{N}}<0$ there is a net energetic stabilization of U, which should have a destabilizing effect on the protein.

Figures 4A and 4B show the ϕ_c-dependence of $E_{\text{cp}}^{\text{U}}$ and $E_{\text{cp}}^{\text{N}}$ for two different values of ϵ_att. For hydrophobic crowders and ϵ_att = 1.0, $E_{\text{cp}}^{\text{N}}$ (open circles) exhibits a slight positive curve and $E_{\text{cp}}^{\text{U}}$ a slight negative curve. These trends are in line with an expanded U with some nonpolar amino acids available for favorable interactions with the crowders, and an N with a hydrophobic core well shielded from the crowders. For ϵ_att = 1.5, the trend is similar but $E_{\text{cp}}^{\text{N}}$ now exhibits a slight negative curve, indicating that, at this strength of attractions, N becomes slightly distorted to accommodate favorable contacts with the hydrophobic crowders. However, these crowder-proteins interactions are rather limited. At the highest packing fraction, ϕ_c = 0.40, $E_{\text{cp}}^{\text{N}}\approx -1.5$, corresponding to ≈ 1–2 fully formed contacts between an nonpolar amino acid and a crowder. In comparison, for the nonspecific case, interactions between crowders and protein are much more prevalent for both U and N. $E_{\text{cp}}^{\text{U}}$ is a sharply decreasing function of ϕ_c for both ϵ_att = 1.0 and 1.5, which will strongly stabilize U at high ϕ_c. However, the decrease in $E_{\text{cp}}^{\text{U}}$ is nearly fully compensated by a decrease in $E_{\text{cp}}^{\text{N}}$. Comparing the two types of crowders, we find that the net energetic effect, $\Delta E_{\text{cp}}=E_{\text{cp}}^{\text{U}}-E_{\text{cp}}^{\text{N}}$, which is what drives destabilization, is more negative for the hydrophobic crowders than for the nonspecific crowders, as can be seen in Figs. 4C and 4D.

Hence, an interesting difference between the two crowder types is that, while their net effects on stability are similar, they differ substantially in the degree of association with the protein. While hydrophobic crowders associate weakly with the protein, even in U, the protein-crowder association in the case of nonspecific interactions is greater by approximately an order of magnitude, as quantified by the magnitude of the interaction energies in the native state, $E_{\text{cp}}^{\text{N}}$ (see Figs. 4A and 4B).

Balancing destabilizing soft interactions and stabilizing steric repulsions

The α₃₅ protein is net stabilized by the presence of hydrophobic crowders at ϵ_att = 1.0 but net destabilized at ϵ_att = 1.5. This suggests a critical attraction strength at which there is no net change in stability. Indeed, for ϵ_att = 1.25, the relative change in T_m is at most one percent and the change in folding free energy, ΔΔF, is very close to zero, as shown in Figs. 5A and 5B. Interestingly, this balance between stabilizing and destabilizing effects is largely independent of ϕ_c. The crowding effect on stability at any given ϵ_att is, however, generally dependent on T, as shown in Fig. 5C for ΔΔF. There is a T-dependence also for excluded volume crowders, which are strictly stabilizing (ΔΔF < 0). Such a T-dependence is possible for excluded volume crowders despite their lack of energetically favorable energetic interactions with the protein, because the size of U is not constant but vary with T (Bazmi & Wallin, 2022). For hydrophobic crowders, due to the energetic stabilization of U relative to N, at low enough temperatures, the net effect crosses over from stabilizing to destabilization, i.e., ΔΔF changes from negative to positive. This crossover temperature (Zhou, 2013), T_c, increases with the contact strength ϵ_att, as shown in Fig. 5D. Our findings for α₃₅ are similar to previous studies on other proteins, such as ubiquitin, which showed a crossover temperature in the presence of either synthetic polymers or protein crowders (Wang et al., 2012b; Zhou, 2013).

Polypeptide crowders

We turn now to the folding of α₃₅ in the presence of polypeptide chains. We use the same model for the peptides as for α₃₅, with the restriction that energetically favorable interactions (i.e., hydrophobic and hydrogen bond interactions) between crowding peptides are turned off. This avoids the formation of oligomeric peptide structures, which would complicate the analysis. We consider two different 5-amino acid sequences: ppppp (peptide 1) and pphpp (peptide 2). Both peptides can thus interact with the α₃₅ chain through backbone-backbone hydrogen bonding and peptide 2, due to its central h amino acid, can additionally interact with α₃₅ through hydrophobic attractions.

Figures 6A and 6B show the T-dependence of the free energy of folding, ΔF, for different numbers N_pep of peptide 1 or peptide 2 chains added to the system. For N_pep ≤ 42, peptide 1 induces a weak stabilization at high T and weak destabilization at low T. This behavior is consistent with a competition between entropy-driven stabilization due to volume exclusion by the peptides and energy-driven destabilization due to inter-chain hydrogen bonding. Interestingly, there is a rather broad temperature range around $T_{\text{m}}^0$ (k_BT ≈ 0.45–0.60) with no detectable change in ΔF (see Figs. 6A and 6D). The apparent lack of crowding effects in this temperature range can arise either because (i) the crowding effects are overall weak at the studied peptide concentrations or (ii) the two opposing crowding effects, peptide excluded volume and peptide-protein attractions, are equal in magnitude and therefore cancel out. To determine which scenario holds, we performed additional simulations with peptides that were only allowed to interact with other chains via repulsive interactions. These “excluded volume peptides” significantly stabilize α₃₅ across all Ts (see Figs. 6C and 6D), which means that scenario (ii) above holds. Hence, in a relatively broad range around the midpoint temperature of the protein, the excluded volume effect due to the peptides is almost perfectly counteracted by soft interactions in the form of hydrogen bonds.

Peptide 2 provides a stronger destabilizing effect on α35 than peptide 1 at a given concentration, as seen in Figs. 6B and 6D. We therefore study up to N_pep = 28 peptide 2 chains. At T ≈ T⁻, ΔF monotonically increases with the number of added peptide 2 chains in contrast to the flat ΔF exhibited by peptide 1. At very low T and N_pep = 28, peptide 2 even leads α₃₅ to become net unstable (ΔF > 0), although the error bars are larger at the lowest studied Ts. A similar behavior is seen for peptide 1 and N_pep = 56. Structurally, low-energy conformations obtained for the peptide 1 and peptide 2 systems exhibit non-native features, including β-sheet structure. An example of a low-energy α₃₅ structure in the precence of peptide 2 crowders is given in Fig. 7A. The destabilization at low T is, at least in part, driven by energetically favorable crowder-peptide interactions that also causes non-native structures. By contrast, low-energy conformations remain native-like in the case of excluded volume peptides (see Fig. 7B). It is instructive to get a rough idea of the volume fraction ϕ_c occupied by the peptides in our systems. Assuming each atom in our C_β-model occupies volume according to its van der Waals radius, we obtain ϕ_c ≈ 0.07 for N_pep = 28. Alternatively, if we assume amino acids are spheres with radius 3.8 Å (typical C_α- C_α distances of peptide bonds) we obtain a slightly higher value, ϕ_c ≈ 0.12. What makes peptide 2 more destabilizing than peptide 1? Because the two peptides are geometrically identical, any difference must derive from soft interactions. Figures 8A and 8B show the peptide concentration dependence of the quantity $\Delta E_{\text{cp}}=E_{\text{cp}}^{\text{U}}-E_{\text{cp}}^{\text{N}}$ (Fig. 4), determined separately for peptide (crowder)-protein hydrogen bond ($\Delta E_{\text{cp}}^{\text{hbond}}$) and hydrophobic ($\Delta E_{\text{cp}}^{\text{hp}}$) energies. For peptide 1, $\Delta E_{\text{cp}}^{\text{hp}}=0$, since this peptide lacks a hydrophobic amino acid. For peptide 2, $\Delta E_{\text{cp}}^{\text{hp}}<0$ for all peptide-protein systems. This means that, vis-à-vis hydrophobic interactions, peptide 2 interacts more favorably with U than with N, causing α₃₅ destabilization. The formation of hydrogen bonds between peptides and the protein are also destabilizing because generally $\Delta E_{\text{cp}}^{\text{hbond}}<0$. Interestingly, $\Delta E_{\text{cp}}^{\text{hbond}}$ is more negative for peptide 2 than for peptide 1, such that the presence of a hydrophobic amino acid on peptide 2 also enhances the formation of hydrogen bonds with the protein chain. In other words, there is a cooperative effect between hydrophobic and hydrogen bond interactions in peptide 2, which further enhances the destabilizing effect of this peptide relative to peptide 1 which has only polar amino acids.