D614G substitution at the hinge region enhances the stability of trimeric SARS-CoV-2 spike protein

DOI: https://doi.org/10.21203/rs.3.rs-154775/v1

Abstract

Background

Spike protein is a key player in the SARS-CoV-2 infection by mediating primary contact between the virus and host cell surface. In the current COVID-19 pandemic, a variant of SARS-CoV-2 having D614G substitution in the spike protein has become dominant world-wide. Initial characterization of the virus shows that the G614 variant is more infectious and has higher fitness than the ancestral (D614) variant. In this study, we analyzed the significance of the D614G substitution on the protein flexibility, inter-residue interaction energies and thermostability of the spike protein trimer.

Results

Using Gaussian network model-based normal mode analysis, we demonstrate that D614G substitution occurs at hinge region that facilitates domain-domain motions between receptor binding domain and S2 region of the spike protein. Further, in-silico mutagenesis and inter-residue energy calculations reveal that contacts involving D614 are energetically frustrated whereas contacts involving G614 are energetically favourable implying the substitution strengthens intra- as well as inter-protomers association. Upon glycine substitution, free energy difference (ΔΔG) is -2.6 kcal/mol for closed and − 2.0 kcal/mol for 1-RBD up conformation i.e., thermodynamic stability has increased. When we perform reverse mutation in the structures of spike protein having G614 substitution, we observe that the free energy difference is 6.6 kcal/mol and 6.3 kcal/mol for closed and 1-RBD up conformations respectively indicating lowered thermodynamic stability. Together, these observations suggest that D614G substitution could modulate the flexibility of spike protein and confer enhanced thermodynamic stability.

Conclusion

Our results on protein flexibility and energetic basis of enhanced stability hint that G614 likely increases the availability of functional form of spike trimer thereby associated to increased infectivity.

Background

According to epidemiological surveillance of the disastrous COVID-19 pandemic, the causing agent severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) virus harbours mutations and is associated with geographical-specific etiological effects [1, 2]. Currently, three major variants of SARS-CoV-2 have been identified namely D614G in the spike protein, G251V in the non-structural protein 3 (NS3) and L84S in the ORF8 protein [3]. The emergence of other variants in recent times from England and South Africa with mutations (e.g., N501Y) in the spike protein are also deemed to be of concern [4–6]. In this article, we focused on the spike protein with D614G substitution. Spike protein of SARS-CoV-2 is a 1273aa long transmembrane glycoprotein and comprises of three modules namely a large ectodomain that protrudes from the surface, a single-pass transmembrane anchor and a short intracellular tail. The ectodomain has S1 and S2 regions responsible for host cell binding and viral-host membrane fusion, respectively. At the junction of S1 and S2 regions, S1/2 cleavage site is present and S2’ cleavage site is located in the S2 region. S1 region comprises of N-terminal domain (NTD, 27–305), receptor binding domain (RBD, 331–528), C-terminal domain 2 (CTD2, 529–590) and C-terminal domain 3 (CTD3, 591–685). Depending on the orientation of RBD in the S1 region, protomers in the functional form of spike protein trimer adopts ‘open’ or ‘closed’ conformation [7–9].

Upon open conformation, RBD exposes ACE2 receptor binding region and interacts with a peptidase domain of angiotensin-converting enzyme 2 (ACE2) receptor. This primary step clasps the virus on to the host surface [10, 11]. Studies on SARS-CoV have shown that subsequent proteolysis at S1/S2 cleavage site sheds S1 region from the spike protein and cleavage at S2’ site near fusion peptide causes a large conformational change in the S2 region. Such conformation change leads to an insertion of fusion peptide to host membrane and the formation of six-helix bundle. At this state, spike protein bridges viral envelope and host membrane. Hairpin-like bend in S2 region brings both membranes to close proximity for fusion following which genetic material gets released into the cytoplasm of the human cell [12]. It is also noted that due to multi-basic nature of S1/2 cleavage site, the SARS-CoV-2 spike protein can be preactivated by furin enzyme during viral packaging [13]. In contrast to SARS-CoV infection, this process reduces the dependence of SARS-CoV-2 on target cell proteases for the succeeding infection.

Therefore, mutations in the spike protein that influence the initial step for viral infection are associated with altered virus transmissibility and pathogenicity [1, 14]. The ancestral spike protein with aspartate at 614th position (SD614) has been asynchronously superseded by glycine substitution (SG614) world-wide. The dominant SG614 variant is shown to have higher infectivity and more efficient replication than the SD614 variant [15–17]. Concomitantly, other studies report that glycine substitution disrupts a salt bridge interaction between aspartate at 614th position and threonine at 859th position of an adjacent protomer and may contribute to higher frequency of open conformation than the SD614 variant [18]. A recent cryo-EM study further reveals that the glycine substitution prevents premature shedding of S1 region while promotes protease cleavage at S1/S2 site [19, 20]. Along the similar line, we demonstrate that the position 614 is present at the hinge region that facilitates domain-domain motions. The substitution of glycine potentially influence the flexibility of hinge region and consequent in the altered RBD motions as observed experimentally [18, 20]. Our calculations of local interaction energies and the free energy difference between SD614 and SG614 variants of the spike protein indicate that glycine creates energetically favourable local environment. As a result, it strengthens the association of S1 and S2 regions of the same as well as adjacent protomer(s) and enhances the overall stability of the spike protein trimer.

Results And Discussion

D614G is located at the hinge region of the spike protein

The conformational transition between closed and open states of spike protein is mediated by large scale motions of RBD. Comparison of closed (PDB code: 6VXX) and 1-RBD up conformation (PDB code: 6VYB) [21] of spike trimer structures shows that apart from RBD that undergoes large movements, NTD and CTD2 domains have considerable Cα deviations with RMSD values of 1.6Å and 3.4Å, respectively (Fig. 1A). However, CTD3 that clasps the position 614 in the S1 region superposes well (RMSD 0.6Å), hinting that it likely facilitates domain-domain motions. To understand the association of 614th position to the domain motions, we performed GNM-NMA for spike protein trimer. In the Gaussian modes, hinges have the characteristics of crossover of displacement profiles from negative to positive or vice versa. The slowest Gaussian modes of protomers in spike protein trimer representing closed and 1-RBD up conformations show that displacement profiles of CTD3 are around zero, indicating the domain is steady while the tethered NTD and RBD exhibit notable displacements (Fig. 1B). Consistently in both conformations, displacement profiles crossover around Lys310-Phe318 and Gly593-Val618 in the S1, indicating they serve as hinges (vertical cyan bars in Fig. 1B). The hinge-1 (Lys-310-Phe318) mediates NTD-RBD motions while the hinge-2 (Gly593-Val618) mediates domain motions between RBD and S2 region (Fig. 1B). Notably, hinge-2 harbours the position of dominant D614G substitution. Given that glycine is a highly flexible residue, the replacement of aspartate to glycine potentially influences the flexibility of hinge-2, which aids the RBD-S2 motion essential for the conformational transition.

D614G substitution strengthens intra- and inter-protomer interactions

Further, we probed the effect of glycine substitution on the energetics of local inter-residue contacts. This can be quantified as frustration of a residue or inter-residue contacts. We calculated frustration index of residues and inter-residue contacts for two variants of spike protein viz. aspartate or glycine at the 614th position. The result shows that frustration index of aspartate in the spike protein (SD614) is -1.25, -1.25 and -1.30 for three protomers in the closed conformation (red lines in Fig. 2A). The frustration index of aspartate in the 1-RBD up conformation is -1.24, -1.31 and -1.28 for three protomers (red lines in Fig. 2B). Hence, in both the conformations aspartate is highly frustrated. Conversely, in glycine variant (SG614), the residue is neutrally frustrated with frustration index of -0.48, -0.42 and -0.46 for protomers in the closed conformation and -0.50, -0.35 and -0.37 for protomers in the 1-RBD up conformation (blue lines in Fig. 2). Recently, structures of spike protein with D614G substitution were released in PDB (November 4, 2020) [20]. When we analyze them, we found that glycine has frustration index of -0.55, -0.19 and -0.88 for protomers in closed conformation and -0.75, -0.39 and -0.31 for protomers in 1-RBD up conformation (PDB codes: 7KDK and 7KDL). These result imply that residue frustration at the 614th position has become neutral upon glycine substitution.

In the spike protein of both conformations (SD614), aspartate is involved in intra-protomer contacts (with residues Ser591, Gly593 and Gly594) as well as in inter-protomer contacts (with Asn616, Arg646, Ser735, Thr859 and Pro862) through direct, long-range electrostatic or water-mediated interactions. The mutational frustration index indicates that all the 8 contacts are highly frustrated (Fig. 3, top panel). However, in the closed conformation of SG614 variant, glycine contacts with Leu611 and Cys649 of the same protomer and Pro862 of the adjacent protomer. Except inter-protomer contact through Pro862, the other intra-protomer contacts are minimally frustrated (Fig. 3, top left panel). Likewise, in the 1-RBD up conformation of SG614 variant, glycine has the same contact pattern, as observed in the closed conformation in addition to contact with Asn616. Of these four contacts, two contacts involving Leu611 and Cys649 are minimally frustrated (Fig. 3, top right panel). Similar observation is seen for mutational indices of contacts involving glycine at 614th position in the cryo-EM structures of D614G variant of spike protein (PDB codes: 7KDK and 7KDL, Supplementary Table S1A). Overall, the number of highly frustrated contacts are reduced upon aspartate to glycine substitution.

Next, we calculated configurational frustration index that indicates how favourable the native contact between two residues relative to other possible contacts those two residues can have. Results show that in the closed conformation, aspartate (SD614) has one minimally frustrated contact with Arg646 (Fig. 3, left bottom panel). Whereas, glycine (SG614) has six minimally frustrated contacts with residues Ser591, Gly593, Asn616, Thr645 and Arg646 of the same protomer and Thr859 of a preceding protomer in the clock-wise direction (Fig. 3, left bottom panel). Similar trend is seen for 1-RBD up conformation in which aspartate (SD614) has a highly frustrated contact with Gly593 while glycine (SG614) has the same contact but minimally frustrated besides three minimally frustrated contacts with other residues (Thr645, Arg646 and Thr859) (Fig. 3, right bottom panel). These observations are common among three protomers present in the spike protein trimer and corroborate with configurational indices calculated from cryo-EM structures of spike protein having D614G substitution (Supplementary Table S1B). Hence, glycine has more favourable contacts than aspartate. Overall, calculations of single residue, mutational and configurational frustrations reveal that glycine substitution modifies local interaction energy in the favourable direction.

D614G substitution increases the thermodynamic stability of spike protein trimer

If the reduction of frustration in the local interaction energies is significant upon aspartate to glycine substitution, it can have an influence on the thermodynamic stability of the spike protein trimer. To examine this, we calculated difference in the total free energy of trimer between SD614 and SG614 variants using FoldX package [22]. Results show that the free energy difference (ΔΔG) is -2.6 kcal/mol for the closed conformation and -2.0 kcal/mol for the 1-RBD up conformation. To affirm this, we performed reverse mutation in the cryo-EM structures of spike protein with D614G substitution (PDB codes: 7KDK and 7KDL). Result shows that the energy differences are positive values i.e., 6.6 kcal/mol and 6.3 kcal/mol for closed and 1-RBD up conformations respectively, indicating destabilizing effect of reverse mutation. These energy differences are higher than the reasonable threshold value of ±0.5 kcal/mol [23]. Hence, the stabilizing effect of glycine substitution in local environment significantly increases the overall stability of spike protein trimer. Together, these results imply that enhanced stability of SG614 may increase the availability of functional form of spike protein trimer and consequent in higher infectivity compared to the SD614 as observed in the recent experimental studies [15,16,19,24].

Conclusions

The increasing severity in public health and economic crisis builds urgency to develop therapeutic intervention against COVID-19 infection at the earliest. Currently, the dominance of D614G substitution in the SARS-CoV-2 spike protein, that is being intensively studied across the globe for COVID-19 prophylaxis and treatment, invites special attention. This study demonstrates using normal mode analysis and in-silico approaches that D614G substitution occurs at the hinge region and potentially influences the conformational transition essential for human ACE2 receptor binding. Glycine changes the local environment from energetically frustrated into favourable for contacts present within as well as between protomer(s). Consequently, the free energy of SG614 is lower than that of SD614, indicating the local changes in the interaction energies at the 614th position in each protomer have a significant effect on the overall thermodynamic stability of the spike protein trimer. This finding bestows to our knowledge on the mechanism of increased transmissibility of SG614.

Methods

Gaussian network model-based normal mode analysis

As spike protein undergoes transition between closed and open conformation involving large scale motions of RBD, it requires assisting regions referred to as ‘hinges’ [25]. To identify regions that precisely act as hinges for domain-domain motions in the S1 region, we performed Gaussian network model (GNM)-based normal mode analysis (NMA) for the spike protein trimer, GNM-NMA is a robust method to accurately predict hinges [26,27]. Based on the approximation that interatomic distances follow Gaussian distribution about their equilibrium values, this method constructs Cα atom-level elastic network model and connects nodes that fall within 7Å distance by virtual spring. With the uniform spring constant for all the connections, Kirchhoff adjacency matrix representing residue proximity is derived. The elements of inverse of the matrix provide information about normal modes of Gaussian dynamics experienced by the spike protein. GNM-NMA calculations were performed for closed (PDB code: 6VXX) and 1-RBD up conformations (PDB code: 6VYB) to increase the confidence of identified hinges. The slowest normal mode was considered for analyzing Gaussian dynamics of the spike protein trimer.

3-D structural model for D614G variant of the spike protein trimer

We generated an in silico model for the D614G variant of spike protein trimer using structure editing tool in UCSF chimera with default parameters [28]. Sidechains were optimized using SCWRL 4.0 program [29]. Two D614G variant models were generated corresponding to closed and 1-RBD up conformations of the spike protein trimer based on the reference cryo-EM structures available in the protein data bank (or PDB) entries 6VXX and 6VYB, respectively. In addition, structures of spike protein with D614G substitution have been recently released [20] and included in this analysis.

Calculation of frustration in the local interaction energy

The effect of D614G substitution on local interaction energies was examined using Frustratometer algorithm [30]. The underlying principle of the algorithm is that a native protein comprises several conflicting residue contacts leading to local frustration. Here, residue contacts are categorized as short range (<6.5Å), long range (6.5-9.5Å) and water-mediated (long range and exposed to solvent) based on Cβ distance. To examine frustrated contacts present in a protein, the algorithm systematically substitutes the residue type or alters chemical configuration of each interacting pair (including water-mediated interactions) and generates approximately 1000 structural decoys for a given contact (elaborated in Ferreiro et al. 2007). The extent of changes in the total interaction energy between native and structural decoys according to associative memory, water-mediated, structure and energy model (AWSEM and implemented as a molecular dynamics algorithm, AWSEM-MD) decides whether the frustration of a concerned contact is minimal, neutral, or high. When the native energy is in the lower end of the energy distribution of structural decoys, the contact is stabilizing and minimally frustrated (favourable), whereas when the native energy falls in between the energy distribution of structural decoys, the contact is neutrally frustrated. Native energy at higher end of the energy distribution of structural decoy indicates that the contact is destabilizing and highly frustrated (unfavourable). Often, highly frustrated contacts signify functional constraints such as substrate binding, allosteric transitions, binding interfaces and conformational dynamics [31,32].

Frustration of the contact is represented as frustration index, a Z-score of interaction energy of native contact with respect to the interaction energy distribution of structural decoys generated for that specific contact. Frustration index below -1 indicates that the interacting pair is highly frustrated while the index between -1 to 0.78 or above 0.78 indicating the interacting pair is neutrally or minimally frustrated, respectively. Depending upon the nature of alteration, frustration is referred as ‘mutational’, ‘configurational’ or ‘single-residue level’. In the mutational frustration, residue type is replaced by other residue types while in the configurational frustration, all possible interaction types between the native residue pairs were sampled through altering residue configuration. In the case of single-residue level frustration, only a single residue is considered. The structural decoy set comprises of randomized residue type at that specific site and frustration index is calculated by evaluating changes in the protein energy upon altering the type of residue. In these three categories of frustration indices, only the concerned site/interaction is altered and the rest of the structure is maintained as native. Here, we analyzed all categories of frustration indices for two variants of spike protein (SD614 and SG614) in the functional trimeric form.

Calculation of thermodynamic stability of SD614 and SG614 variants

To study the effect of D614G variation on the thermodynamic stability of the spike protein trimer, we calculated free energy changes upon aspartate to glycine substitution using buildmodel function in FoldX [22]. Five iterations of free energy calculations were carried out to obtain converged results [33]. Inferences of the results were derived from closed and 1-RBD up conformations of the spike protein trimer.

Abbreviations

ACE2, angiotensin-converting enzyme 2

AWSEM, associative memory, water-mediated, structure and energy model

GNM, Gaussian network model

NMA, normal mode analysis

PDB, protein data bank

RBD, receptor binding domain

SARS-CoV-2, severe acute respiratory syndrome coronavirus-2

SD614, spike protein trimer with aspartate at 614th position

SG614, spike protein trimer with glycine at 614th position

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

All authors have given their consent to publish the findings in this paper.

Availability of data and materials

Supporting data are enclosed as Supplementary Table S1. Data generated or analyzed during this study are included in this published article and its supplementary information file.

Competing interests

The authors declare they have no competing interests.

Funding

This research is supported by the Department of Biotechnology, Government of India through the IISc-DBT partnership program and FIST program sponsored by the Department of Science and Technology, Government of India and also the Council of Scientific and Industrial Research, India. Support from UGC, India – Centre for Advanced Studies and Ministry of Human Resource Development, India is gratefully acknowledged. NS is a J. C. Bose National Fellow.

Author contributions

AY and NS conceived the project and NS supervised the work. AY and DS performed the analysis. AY wrote the draft version of the manuscript and finalized with the help of DS and NS. All authors read and approved of the manuscript.

Author information

Affiliations

Lab 103, Molecular Biophysics Unit, Indian Institute of Science, Bangalore, Karnataka, 560012, India

Arangasamy Yazhini, Das Swayam Prakash Sidhanta and Narayanaswamy Srinivasan

Corresponding authors

Correspondence to Narayanaswamy Srinivasan

References

  1. Brufsky A. Distinct viral clades of SARS-CoV-2: Implications for modeling of viral spread. J. Med. Virol. 2020. p. 1386–90.
  2. Mercatelli D, Giorgi FM. Geographic and Genomic Distribution of SARS-CoV-2 Mutations. Front Microbiol. 2020;11.
  3. Forster P, Forster L, Renfrew C, Forster M. Phylogenetic network analysis of SARS-CoV-2 genomes. Proc Natl Acad Sci USA. 2020;117:9241–3.
  4. Passos JCSGA. The high infectivity of SARS-CoV-2 B.1.1.7 is associated with increased interaction force between Spike-ACE2 caused by the viral N501Y mutation. bioRxiv. 2021;501:1–9.
  5. Luan B, Wang H, Huynh T. Molecular Mechanism of the N501Y Mutation for Enhanced Binding between SARS-CoV-2’s Spike Protein and Human ACE2 Receptor. bioRxiv. 2021;2021.01.04.425316.
  6. Tegally H, Wilkinson E, Giovanetti M, Iranzadeh A, Fonseca V, Giandhari J, et al. Emergence and rapid spread of a new severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2) lineage with multiple spike mutations in South Africa. medRxiv. 2020;10:2020.12.21.20248640.
  7. Walls AC, Park YJ, Tortorici MA, Wall A, McGuire AT, Veesler D. Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein. Cell. 2020;181:281–92.
  8. Wrapp D, Wang N, Corbett KS, Goldsmith JA, Hsieh C-L, Abiona O, et al. Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation. Science. 2020;367:1260–3.
  9. Cai Y, Zhang J, Xiao T, Peng H, Sterling SM, Walsh RM, et al. Distinct conformational states of SARS-CoV-2 spike protein. Science. American Association for the Advancement of Science; 2020;eabd4251.
  10. Lan J, Ge J, Yu J, Shan S, Zhou H, Fan S, et al. Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor. Nature. 2020;581:215–20.
  11. Yan R, Zhang Y, Li Y, Xia L, Guo Y, Zhou Q. Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2. Science. 2020;367:1444–8.
  12. Belouzard S, Millet JK, Licitra BN, Whittaker GR. Mechanisms of coronavirus cell entry mediated by the viral spike protein. Viruses. 2012. p. 1011–33.
  13. Shang J, Wan Y, Luo C, Ye G, Geng Q, Auerbach A, et al. Cell entry mechanisms of SARS-CoV-2. Proc Natl Acad Sci USA. 2020;117:11727–11734.
  14. Li Q, Wu J, Nie J, Zhang L, Hao H, Liu S, et al. The Impact of Mutations in SARS-CoV-2 Spike on Viral Infectivity and Antigenicity. Cell. 2020;182:1284-1294.e9.
  15. Korber B, Fischer WM, Gnanakaran S, Yoon H, Theiler J, Abfalterer W, et al. Tracking Changes in SARS-CoV-2 Spike: Evidence that D614G Increases Infectivity of the COVID-19 Virus. Cell. 2020;182:812-827.e19.
  16. Plante JA, Liu Y, Liu J, Xia H, Johnson BA, Lokugamage KG, et al. Spike mutation D614G alters SARS-CoV-2 fitness. Nature. 2020;1–9.
  17. Hou YJ, Chiba S, Halfmann P, Ehre C, Kuroda M, Dinnon KH, et al. SARS-CoV-2 D614G variant exhibits efficient replication ex vivo and transmission in vivo. Science. 2020;370:1464–8.
  18. Yurkovetskiy L, Wang X, Pascal KE, Tomkins-Tinch C, Nyalile TP, Wang Y, et al. Structural and Functional Analysis of the D614G SARS-CoV-2 Spike Protein Variant. Cell. 2020;183:739-751.e8.
  19. Zhang L, Jackson CB, Mou H, Ojha A, Peng H, Quinlan BD, et al. SARS-CoV-2 spike-protein D614G mutation increases virion spike density and infectivity. Nat Commun. 2020;11:1–9.
  20. Gobeil SM-C, Janowska K, McDowell S, Mansouri K, Parks R, Manne K, et al. D614G mutation alters SARS-CoV-2 spike conformation and enhances protease cleavage at the S1/S2 junction. Cell Rep. 2020;108630.
  21. Schrödinger. The PyMol (TM) Moelcular Graphics System, Version 1.7.2.1 Schrödinger, LLC. 2015;
  22. Schymkowitz J, Borg J, Stricher F, Nys R, Rousseau F, Serrano L. The FoldX web server: An online force field. Nucleic Acids Res. 2005;1:W382-388.
  23. Bromberg Y, Rost B. Correlating protein function and stability through the analysis of single amino acid substitutions. BMC Bioinformatics. 2009;10.
  24. Zhang J, Cai Y, Xiao T, Lu J, Peng H, Sterling SM, et al. Structural impact on SARS-CoV-2 spike protein by D614G substitution. bioRxiv. 2020;2020.10.13.337980.
  25. Bahar I, Lezon TR, Bakan A, Shrivastava IH. Normal mode analysis of biomolecular structures: Functional mechanisms of membrane proteins. Chem Rev. 2010;110:1463–97.
  26. Haliloglu T, Bahar I, Erman B. Gaussian Dynamics of Folded Proteins. Physcial Rev Lett. 1997;79:3090.
  27. Bahar I, Atilgan AR, Erman B. Direct evaluation of thermal fluctuations in proteins using a single-parameter harmonic potential. Fold Des. 1997;2:173–81.
  28. Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, et al. UCSF Chimera - A visualization system for exploratory research and analysis. J Comput Chem. 2004;25:1605–12.
  29. Krivov GG, Shapovalov M V., Dunbrack RL. Improved prediction of protein side-chain conformations with SCWRL4. Proteins Struct Funct Bioinforma. 2009;77:778–95.
  30. Parra RG, Schafer NP, Radusky LG, Tsai MY, Guzovsky AB, Wolynes PG, et al. Protein Frustratometer 2: a tool to localize energetic frustration in protein molecules, now with electrostatics. Nucleic Acids Res. 2016;44:W356–60.
  31. Ferreiro DU, Hegler JA, Komives EA, Wolynes PG. Localizing frustration in native proteins and protein assemblies. Proc Natl Acad Sci USA. 2007;104:19819–24.
  32. Ferreiro DU, Komives EA, Wolynes PG. Frustration in biomolecules. Q. Rev. Biophys. 2014. p. 285–363.
  33. Tokuriki N, Stricher F, Schymkowitz J, Serrano L, Tawfik DS. The Stability Effects of Protein Mutations Appear to be Universally Distributed. J Mol Biol. 2007;369:1318–32.