Protein tertiary structure and the myoglobin phase diagram.

We develop an effective theory approach to investigate the phase properties of globular proteins. Instead of interactions between individual atoms or localized interaction centers, the approach builds directly on the tertiary structure of a protein. As an example we construct the phase diagram of (apo)myoglobin with temperature (T) and acidity (pH) as the thermodynamical variables. We describe how myoglobin unfolds from the native folded state to a random coil when temperature and acidity increase. We confirm the presence of two molten globule folding intermediates, and we predict an abrupt transition between the two when acidity changes. When temperature further increases we find that the abrupt transition line between the two molten globule states terminates at a tricritical point, where the helical structures fade away. Our results also suggest that the ligand entry and exit is driven by large scale collective motions that destabilize the myoglobin F-helix.

Introduction

In the description of a complex system such as a protein, it is often impractical, if not impossible, to accurately model physical phenomena with a “fundamental” level precision. For example, practical chemistry, including even precision first-principles quantum chemistry, is never concerned with the detailed structure of the atomic nucleus. Instead, one focuses on a few key variables and constructs an effective theory for those. In many circumstances and in particular when the system admits either symmetries or a separation of scales, the reduced set of variables can then be treated on its own right.

In the case of proteins, the great success of structural classification schemes such as SCOP[1] and CATH[2] and many others[3], demonstrates that folded proteins are built in a modular fashion, and from a relatively small number of components that are made of several amino acids. Here we exploit this modularity of protein structures to construct the phase diagram of globular proteins, with temperature (T) and acidity (pH) as the thermodynamical variables. The methodology that we develop is very general and applicable to most globular proteins[4], even though we here develop it using myoglobin (Mg)[5-7] as a concrete example: Our approach is based on the Landau-Ginsburg-Wilson (LGW) paradigm[8,9] which is a systematic way to construct effective theories that model the properties of different phases in a material system, and transitions between them. Instead of individual atoms or other highly localized interaction centers and their mutual interactions, our effective theory description employs the entire tertiary structure of a protein as the fundamental constituent: We describe the protein backbone as a single multi-soliton[10,11]. This multi-soliton acts as an attractor in the energy landscape, it is a minimum free energy state towards which other conformations become funneled. Indeed, a multi-soliton solution to a non-linear difference (differential) equation is the paradigm structural self-organizing principle in many physical scenarios. Here it emerges as a stable solution to a variational equation that we obtain from the LGW free energy, and it governs the mutual interactions between the individual solitons that model the super-secondary structures (helix-loop-helix etc.) of the protein.

The advantage of the LGW formalism in combination with the soliton-concept is computational efficiency, over any other approach to protein dynamics that we are aware of; the method enables us to perform numerical simulations and analyses with very high efficiency. Similar approaches have proven highly successful in many complex scenarios with extended filament-like objects, from the description of interacting superconducting vortex lines to complex knotted fluxtubes[12,13]. Indeed, the evaluation of a T−pH phase diagram of any protein using e.g. molecular dynamics would be unthinkable, with presently available computers.

Myoglobin is a stable, relatively simple globular protein that is the paradigm example in protein folding and unfolding studies[14-30]. It plays an important role in biological processes such as electron transfer, oxygen delivery, catalysis and signaling. In particular, myoglobin can bind small non-polar ligands such as O2, CO, and NO in its interior, where they become attached to the iron atom of the heme. The native folded state (N) of myoglobin is very compact, and supports eight α-helices (labelled A to H). Since there is no apparent static channel for the ligands to enter and exit, myoglobin must undergo conformational deformations for the ligands to pass[5-7,31-36]. These deformations are regulated by changes in physiological conditions in particular by variations in temperature and/or acidity. Several experiments have been performed to investigate (un)folding pattern of myoglobin as a function of pH, mostly at room temperature. These experiments reveal that the (un)folding proceeds reversibly and sequentially, according to a four-state scheme: At low pH values near or below around pH~2 the structure resembles a random coil (U). In the vicinity of the regime two folding intermediates I and I can be found[18,19], both akin a molten globule[37,38] with a structure that changes with varying pH. When pH reaches values that are above 4–4.5 (apo)Mg starts entering its native folded state. Overall, the transitions seems to proceed according to the scheme N ↔ I ↔ I ↔ U as the acidity changes[18]. Variable-T experiments are less common, but the results are quite similar[27-29]: At very low temperatures and near neutral pH, the structure of Mg is in the folded native state N; the F-helix of apoMg is disordered, even at relatively low temperatures. Below T~340 K and close to neutral pH the structure remains in the vicinity of its native state. When T increases further Mg becomes a molten globule, and when T becomes even higher the helical content starts to decrease: The first to unfold is helix F followed by helices B,C,D and E. Then the helices A, G and H loose their stability. At very high temperatures the structure resembles a random coil.

Methods

Our effective theory approach describes the (virtual) Cα protein backbone in terms of the (virtual) bond (θ) and torsion (ϕ) angles. To evaluate these coordinates, we frame the Cα backbone by the mutually orthonormal backbone tangent (t), binormal (b) and normal (n) vectorswhere r (i = 1,..., n) are the Cα coordinates. These vectors are subject to the discrete (Frenet) equation[39]Here T2 and T3 generate three dimensional rotations, with (T) = ε. From (1), (2) we can determine (θ, ϕ) in terms of the Cα coordinates r. Conversely, when θ and ϕ are all known we can reconstruct the r by solving (2) (for details see Supplementary Material) when we assume that the distance between neighboring Cα atoms remains close to the average value ~3.8 Å: A good quality all-atom approximation of the entire heavy atom structure of a protein can always be reconstructed from the knowledge of the (θ, ϕ) coordinates[3,40-43]. In particular, we may employ (θ, ϕ) as the variables in a free energy that models the protein backbone.

Previously, a number of effective energy functions for the Cα backbone have been constructed using the coordinates (θ, ϕ). Familiar examples include the fully flexible chain model and its extensions[44-46], that are widely used in studies of biological macromolecules and other filament-like objects.

Here we introduce a free energy description that is designed to model folded proteins and their properties at the level of the tertiary structures[39,47-51]. The structure of the energy landscape is determined by the following free energy (for details see Supplementary Material)

In the first sum of (3) we recognize the structure of the energy function of the discretized non-linear Schrödinger (DNLS) equation in the Hasimoto representation[10]. The second sum of (3) then extends the DNLS energy function so that it model folded proteins: The first two terms in the second sum are both among the conserved charges in the DNLS hierarchy, they are called the momentum and the helicity respectively. Both of these terms are odd in torsion angles, thus they break parity which makes the backbone (right-handed) chiral. The third term of the second sum is the Proca mass, together with the second term in the first sum it comprises the original Kirchhoff energy of an elastic rod[45]. Finally, the last term is a hard-core Pauli repulsion with a step-wise profile, it ensures that the distance between any two Cα atoms is at least 3.8 Å (for detailed analysis of V(r) see[51]).

Note that there is no need to introduce any long distance contribution to V(r). The long distance interactions are already accounted for by the properties of the solution to the extended DNLS equation: The DNLS equation is the prototype integrable difference equation that supports solitons as its classical solutions. Solitons are the paradigm examples of extended self-organized objects in a physical system[11]. For appropriate parameter values the DNLS free energy (3) models the entire tertiary structure of a given folded protein backbone, as a single stable minimum energy multi-soliton solution to the variational equations

The multi-soliton profile then describes the various super-secondary structures such as helix-loop-helix (regular-loop-regular) as mutually interacting individual solitons[52,53]. Over a single soliton profile the parameter values in (3) are uniform, and since a soliton extends over several amino acids the number of parameters is generically much smaller than the number of amino acids. In the case of a myoglobin, we use the Protein Data Bank (PDB) structure 1ABS (sperm whale)[33] as a decoy to construct the parameters. This structure has been measured at the very low liquid helium temperature value of around 20 Kelvin and as a consequence the thermal B-factors are very small. We identify ten individual DNLS solitons profiles along the 154 residue 1ABS backbone that become combined into a single multi-soliton solution of the DNLS equations (4), (5) with ~0.8 Å Cα root-mean-square-distance (RMSD) precision[47] (see also Supplementary Material).

We construct the T-pH phase diagram by computing the grand canonical ensemble of statistical physics, and evaluate the observables by averaging them over all possible tertiary structures, weighted by the grand canonical distribution with free energy (3):Here is a normalization factor, β is the inverse temperature factor and μN is a chemical potential contribution to be specified. Since the trace extends over all possible tertiary structures, the entropy contribution relates to the number of all possible tertiary structures. We evaluate (6) numerically, using the Glauber algorithm with acceptance ratio determined by the probability distribution[51,54,55]Here is the variation of between consecutive Monte Carlo steps. We note that Glauber algorithm models pure relaxation dynamics, and for simple systems it reproduces Arrhenius law. At the same time it has been found that small proteins fold according to Arrhenius law[56]. We also note that the (inverse of the) Glauber temperature factor β does not coincide with physical temperature factor kT where k is the Boltzmann constant and T is measured in Kelvin’s, instead the relation is determined by renormalisation group techniques[47,57].

We determine the chemical potential contribution in (6), by recalling the Henderson-Hasselbalch equation[58] that relates the concentrations of protonated and non-protonated amino acids to the difference between pH and acid dissociation constant pK. On the other hand, Gibbs free energy is commonly taken to vary with acidity as follows,where is the protonation probability of a particular amino acid. Histidine with pK~6.0 is the only amino acid in the genetic code that has strong reactivity to pH variations in the physiologically important range from pH~8 down to pH~4. For lower pH both glutamic acid (pK~4.2) and aspartic acid (pK~3.9) need to be accounted for. For simplicity, here we only aim to model the phase diagram for pH above ~4 and up to neutral value so that we only need to account to the contribution of the histidines in 1ABS. Thenand to leading order

We recognize in (9) the format of the Glauber transition probability (7). Moreover, since the DNLS hierarchy admits a unique conserved number operator N~θ2 [10] we propose that in the LGW approachwhere the summation extends over the histidines of 1ABS. As a consequence, to leading order in the LGW approximation μ depends linearly on pH. We also note that for 1ABS pH = 9.0. Accordingly we normalize μ = 0 at that value, to ensure that the ensuing multi-soliton profile models the 1ABS backbone.

As order parameters a.k.a. reaction coordinates we use the radius of gyration R and the α-helical content . We compute their (T, μ) dependence numerically from (6). We take a Cα atom to be in an α-helical posture when for (θ, ϕ) both |θ−θ0| ≤ 0.14 (rad) and |ϕ−ϕ0| < 0.3 (rad) where θ0 = 1.55 and ϕ0 = 0.9 are the PDB average values of the α-helical bond and torsion angle. The counts the relative number of residues in α-helical posture as a function of T and pH; most PDB myoglobins have a value 72–78%, and for 1ABS .

We have simulated 5.000 independent heating and cooling (unfolding and folding) trajectories using the Glauber algorithm, obtained by varying the (inverse) temperature factor β; the trajectories are equally distributed between 50 values of μ ∈ [0, 0.05]. Along each trajectory we first increase temperature (i.e. decrease β) at an adiabatically slow rate, so that the system remains very close to a thermal equilibrium for all β. The value of β is also kept at its high temperature value for a large number of simulation steps, for full thermalization. Finally, the system is brought back to the low temperature value, by reversal of the heating procedure: We have been extremely careful to always thermalize the ensemble before we evaluate any observable[51].

Results

In Fig. 1 we compare the μ = 0 temperature dependence of the observable to experimentally measured α-helicity of (horse heart) myoglobin during thermal denaturation. The experimental data is adapted from[28]. We use this Figure to relate the Glauber temperature T to Celsius scale. Accordingly, our simulations cover the range 0 °C–120 °C at each μ value. The Figs 2–5 summarizes our findings:

Figure 1

Comparison between simulated (sperm whale) value that counts the relative number of residues in α-helical posture, and experimentally determined (horse heart) α-helical content during thermal denaturation. Experimental data is adapted from[28].

Figure 2

(Dis)ordering temperatures for the eight helices A-H at μ = 0, in terms of deviations in torsion angles from their α-helical values. Note that helix-F becomes disordered at relatively low temperature.

Figure 5

The phase diagram in terms of the observable that counts the relative number of residues in α-helical posture, at low temperatures T < 40 °C.

Fig. 2 shows the helix (dis)ordering during a heating and cooling simulation cycle, as a function of temperature at μ = 0 and in terms of the average value of torsion angles ϕ; we recall that for an α-helix ϕ ≈ 1 (rad).

We observe that F-helix starts to disorder soon after T = 20 °C and becomes fully randomized slightly above T = 40 °C. The next to disorder are the helices B, C, D and E; this occurs near T = 90 °C. This is followed by disordering of the helices H and A, and the helix G is last to disorder. Slightly above T = 100 °C the entire chain is fully randomized; according to Fig. 1, at these temperature values also reaches its high temperature asymptotic value. All these simulation results are fully in line with experimental observations[24].

Fig. 3 identifies the simulated phase structure on the (T, μ) plane in terms of radius of gyration R. For μ ≈ 0 we confirm the findings of the Figs 1 and 2: The native state (N) is a region with low temperature (T < 30 °C) and very small values μ < 0.003. Beyond this there is a region where the F-helix (dis)orders (F), it extends to around T ≈ 40 °C and to μ-values up to μ ≈ 0.01. When T and μ increase further, we identify a phase that we denote I and identify as a molten globule intermediate; the radius of gyration values are above 20 Å but below 28 Å in this region, depending on values of T and μ.

Figure 3

The R phase diagram: Collapsed native state (N), (dis)ordering of F-helix (F), two molten globules (I and I) and random coil (U) phase are identified. Note the presence of an apparent tricritical point, when the transition line between I and I terminates in U.

We propose that the high sensitivity of the F-helix that we observe, when either temperature or acidity increases from their low values, controls the ligand entry and exit: The F-helix contains the proximal histidine that is connected to the heme. Thus disordering of the F-helix may expose the heme, for ligand transport.

At μ ≈ 0.027 and for values T < 80 °C we observe a rapid transition: The radius of gyration decreases in a jump-like fashion, by around 4–6 Å depending on T. We interpret this to be a transition between two molten globule intermediates so that for μ > 0.027 we have the second molten globule I[18,19].

Most notably, we observe the presence of an apparent tricritical point, in conjuction with the two molten globule states: The transition line between I and I terminates at around T ≈ 80 °C when both molten globules simultaneously enter the random coil phase.

The different regions of the phase diagram in Fig. 3 can be scrutinized using the detailed R values. In Fig. 4 we show how R varies as a function of μ, at T = 0 °C. We observe a clear change in the derivate of R w.r.t. μ at around μ ≈ 0.003 and also around μ ≈ 0.01. These correspond to the transitions between the N and F, and between the F and I regions in the phase diagram on Fig. 2. Note that the R value of the molten globule I is very close to the experimentally reported value R~23.6 Å[16,21,24].

Figure 4

The evolution of R at temperature close to T = 0 °C. The native state, the state with disordered F-helix and the two molten globule states I and I are all identifiable.

We also have a jump-like (discontinuous) transition in R values at around μ = 0.027, between the two molten globules I and I.

Finally, in Fig. 5 we display the values of that counts the relative number of residues in α-helical posture, for T < 40 °C. We observe that there is only a weak dependence on μ, even though we do note a slight change in the overall stability of helix-F even at relatively small μ values. We conclude that at low temperatures the increase of μ appears to have a stronger influence on loops than on helices. In particular, the ligand transport mechanism if indeed associated with instability in helix-F, appears to engage the adjacent loop structures as well.

The results in Fig. 5 are consistent with room temperature CD helicity measurements that report only minor signal variations for pH values above ~4.5[23]: Acidity does not have a strong effect on the hydrogen bonds that stabilize the helical structures. In the Figure we identify the transition between I and I, in terms of a region with (slightly) decreased value of . It is also notable that right prior to the transition, there is region in I with an enhanced values.

Discussion

In summary, we have proposed to model protein thermodynamics directly at the tertiary level of structures, in terms of the multi-soliton solution of the DNLS equation. We have numerically evaluated the ensuing grand canonical partition function at finite temperature and chemical potential, with the latter identified by comparison with the Henderson-Hasselbalch equation. As an example we have constructed the (T, pH) phase diagram of myoglobin. All our results are in a good agreement with experimental observations. In particular, the ordering of helix stabilization and the emergence of two molten globules are qualitatively in full agreement with experimental observations. Furthermore, we observe that the F-helix with its proximal histidine, is the first to loose stability as either temperature or acidity increase from neutral values. This supports that the destabilization of the F-helix region might have a pivotal role for ligand entry and exit. We have also made predictions for future experiments, in particular we have proposed that at high temperatures near T = 80 °C there is an apparent tricritical point where the two molten globules come together with the random coil phase. Our results show that effective theories that model protein structure directly at the tertiary structure level, can provide a viable computational approach to investigate the phase structure of complex globular proteins.

Supplementary Dataset 1