Alpha helices are more evolutionarily robust to environmental perturbations than beta sheets: Bayesian learning and statistical mechanics to protein evolution
Tomoei Takahashi George Chikenji Kei Tokita Yoshiyuki Kabashima


APS/123-QED







takahashi-tomoei@g.ecc.u-tokyo.ac.jp
Institute for Physics of Intelligence, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan.

 

Graduate School of Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, 464-8603, Japan.
 






 
Graduate School of Informatics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, 464-8601, Japan.


Institute for Physics of Intelligence, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan.
Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
Trans-Scale Quantum Science Institute, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

 
 

 










             

§ ABSTRACT 
How typical elements that shape organisms, such as protein secondary structures, have evolved, or how evolutionarily susceptible/resistant they are to environmental changes, are significant issues in evolutionary biology, structural biology, and biophysics. According to Darwinian evolution, natural selection and genetic mutations are the primary drivers of biological evolution. However, the concept of “robustness of the phenotype to environmental perturbations across successive generations," which seems crucial from the perspective of natural selection, has not been formalized or analyzed. In this study, through Bayesian learning and statistical mechanics we formalize the stability of the free energy in the space of amino acid sequences that can design particular protein structure against perturbations of the chemical potential of water surrounding a protein as such robustness. This evolutionary stability is defined as a decreasing function of a quantity analogous to the susceptibility in the statistical mechanics of magnetic bodies specific to the amino acid sequence of a protein. Consequently, in a two-dimensional square lattice protein model composed of 36 residues, we found that as we increase the stability of the free energy against perturbations in environmental conditions, the structural space shows a steep step-like reduction. Furthermore, lattice protein structures with higher stability against perturbations in environmental conditions tend to have a higher proportion of α-helices and a lower proportion of β-sheets. The latter result shows that protein structures rich in α-helices are more robust to environmental perturbations through successive generations than those rich in β-sheets.




                              
Alpha helices are more evolutionarily robust to environmental perturbations
than beta sheets: Bayesian learning and statistical mechanics to protein evolution
    Yoshiyuki Kabashima
    September 9, 2024
==============================================================================================================================================================






§ INTRODUCTION


Understanding whether fundamental elements that shape organisms are evolutionarily robust or prone to change is crucial for addressing why the phenotypes of existing organisms are so limited compared to the physically possible patterns, or for resolving significant issues such as predicting evolution.

Many studies on protein evolution suggest that mutational robustness may have driven the evolution of characteristic protein structures<cit.>. It has been shown that there is a correlation between mutational robustness and the geometric symmetry of protein structures<cit.>, that secondary structures are more robust to mutations than intrinsic disorder structures<cit.>, and that mutational robustness and the structural modularity of proteins (i.e., the proportion of amino acid residues forming secondary structures) contribute to the evolvability of proteins<cit.>. Furthermore, the low algorithmic complexity of gene sequences achieves symmetric protein structures<cit.>. It is also demonstrated that mutational robustness is well-compatible with the functional sensitivity of proteins<cit.>, and there is a correlation between the dynamics of protein structures and their evolvability<cit.>. These various lines of research strongly suggest that mutational robustness (and additionally, the low algorithmic complexity of genes) drives the formation of protein secondary structures and the evolvability of proteins.

In addressing the evolution problem, a statistical mechanics approach to abstract models has also produced significant results. A study using the spin glass model has shown that, in evolved organisms, plasticity due to environmental fluctuations and plasticity due to mutations are strongly correlated, leading to a dimensional reduction in the phenotypic space<cit.>. Studies on gene regulatory network (GRN) models have revealed that GRNs with high fitness exhibit high mutational robustness. This tendency is stronger in GRNs obtained through evolutionary simulations than those obtained through efficient sampling methods for exploring high-fitness GRNs<cit.>. It has also been found that both mutational robustness and developmental robustness drive the evolution of GRNs<cit.>. These studies suggest that biological evolution possesses mechanisms different from typical optimization processes, enabling the selection of GRNs (phenotypes) with high mutational robustness and noise during development. These findings imply the reduction of phenotypic space.

The above results have significantly advanced our understanding of the relationship between mutational robustness, developmental robustness, and evolvability in proteins (or life in general). However, these concepts of robustness pertain solely to the stability of an individual (or phenotype) against various perturbations over its lifetime. Since evolution involves changes in genetic information and associated traits over generations, considering the stability in terms of how well a trait (phenotype) adapts to or is maladapted to its environment, how this influences the genotype that produces the phenotype, and how these influences alter traits in subsequent generations, helps understand evolution.

In this study, we define the free energy of the space of amino acid sequences (i.e., genetic information) that can design a particular protein structure, utilizing the framework of Bayesian learning. We discuss the stability of this free energy against perturbations in environmental conditions surrounding the protein. The stability of the free energy is determined for a randomly generated two-dimensional (2D) lattice protein model, and we elucidate the relationship between the structural features of the protein and the stability of its free energy.

We use lattice proteins for model proteins<cit.>. Random structural patterns that have not evolved do not exist in protein structure databases such as the Protein Data Bank (PDB)<cit.>. Therefore, artificial models like lattice proteins are more suitable for this study. We show the definition of secondary structures in 2D lattice proteins of this study in <ref>. In the analysis of secondary structures of 2D lattice proteins, β-sheets decrease as the designability (the number of amino acid sequences that fold into a given protein structure) increases<cit.>. Lattice proteins are also effective for analyzing the free energy landscape of protein folding<cit.>, the phenomenon of cold denaturation where proteins denature at low temperatures<cit.>, and the impact of amino acid residue mutations on the native structure of proteins<cit.>. Additionally, lattice models are used to analyze the folding energy landscape of RNA<cit.>. Therefore, if one develops a valid theory for protein evolution, it could be said that lattice proteins can reveal qualitatively accurate behaviors in analyzing the environmental robustness of protein secondary structures across successive generations, our objective in this study.



§ MODEL AND METHOD



 §.§ Hamiltonian of lattice HP model with the water chemical potential


The lattice HP model places amino acids at lattice points, representing the protein structure as a self-avoiding walk on the lattice. A self-avoiding walk is a path that does not pass through the same point more than once on a lattice (or graph). The naturally occurring 20 types of amino acids are simplified into two types: hydrophobic (H), which repels water molecules, and polar (P), which attracts water molecules.

The structure of a protein, denoted as R, is represented by the set of coordinates r_i for each amino acid. For a protein with N amino acids, R =  r_1,  r_2, ⋯,  r_N. The state of the i-th amino acid is denoted by σ_i, with σ_i = 1 (for hydrophobic) and 0 (for polar). The lattice HP model typically considers only the attractive interactions between hydrophobic amino acids. However, since proteins also interact with water molecules surrounding them, in this study, the Hamiltonian of a protein with structure R and sequence σ is expressed as follows, with μ representing the chemical potential of water near the protein surface:

    H( R , σ; μ) = -∑_i<jσ_iσ_jΔ( r_i -  r_j) - μ∑_i = 1^N (1 - σ_i),

where Δ( r_i -  r_j) is a “contact function" that equals 1 when the i-th and j-th amino acids are spatially nearest neighbors but not consecutive in the sequence, and 0 otherwise. The second term represents the interaction of polar amino acids with surrounding water molecules. 

This Hamiltonian, proposed for the first time in our previous study<cit.>, includes a hydration term, as shown in Eq. (<ref>), which is crucial for exploring the environmental robustness of protein structures<cit.>. Additionally, Eq. (<ref>) can be seen as a simplified version of the protein energy in water, excluding interactions between water around the protein and bulk water, as shown in studies like<cit.>. It is important to note that μ is not a parameter representing the overall environment within an organism but rather the environment surrounding a specific protein.

The chemical potential of water surrounding a protein, μ, can be considered a parameter representing environmental conditions, as it depends on the state of bulk water (such as pH and pressure) and temperature. The strength of hydrogen bonds between hydrophilic amino acid residues and water molecules is generally influenced by temperature. Consequently, μ can be regarded as an environmental condition specific to the protein structure.



 §.§ Bayesian learning framework

Bayesian learning is a framework for machine learning based on Bayesian statistics, in which one updates the prior probability (prior) of an event to a posterior probability (posterior) in light of observed data. Statistical models derived from Bayesian learning are often interpreted as probabilistic generative models of observational data. We utilize the properties of Bayesian learning to construct a probabilistic generative model for the phenotype of proteins, namely their native structures. The probabilistic generative model of protein structures described below was devised in our previous work as a method for protein design and has achieved a certain level of success in the context of lattice protein models<cit.>. Protein design problem is the inverse problem of protein structure prediction. Protein design is thus determining the sequence of amino acids that will fold into a given protein structure <cit.>.

In our model, we consider the native structure R as the observed variable, the amino acid sequence σ as the latent variable, and the chemical potential of the surrounding water μ as the hyperparameter. In this context, Bayes' theorem is expressed as follows:

    p(σ |  R, μ) = p( R|σ, μ) p(σ|μ)/∑_σp( R|σ, μ) p(σ|μ)

Let β the inverse temperature of the environment, the likelihood function p( R|σ), prior p(σ|μ), and posterior p(σ |  R, μ) are respectively given as

    p(R | σ, μ) =   e^-β H(R , σ ; μ)/Z(σ ; β, μ),
    
         p(σ | μ) =   Z(σ; β_p, μ_p)/Ξ(β_p, μ_p),
    
        p(σ| R, μ) =   e^-β H(R, σ ; μ) /Y( R;β, μ).

The normalization constants of Eqs. (<ref>), (<ref>), and (<ref>) are as follows:

    Z(σ ; β, μ)    = ∑_ R e^-β H(R , σ ; μ),
    Ξ(β_p, μ_p)    = ∑_σ∑_ R e ^ - β_p H(R, σ; μ_p),
    
        Y( R;β, μ)    = ∑_σ e^-β H(R, σ ; μ).

We refer to Eqs. (<ref>), (<ref>), and (<ref>) as the structural partition function, grand partition function, and sequence partition function, respectively. Additionally, the inverse temperature and chemical potential in the prior Eq. (<ref>) may differ from those of protein folding and are thus denoted as β_p and μ_p, respectively.

An important point is that the posterior Eq.(<ref>) does not include the structural partition function Z(σ ; β, μ). This is because Z(σ ; β, μ) requires an exhaustive structural search of the Boltzmann factor, which is infeasible considering the infinite degrees of freedom of protein structures. The sum over sequences ∑_σ is much more manageable than the sum over structures ∑_ R.

The likelihood function p( R | σ, μ) represents the probability of a given structure R occurring for a given amino acid sequence σ. This is the probability that σ folds into structure R. More generally, it is the probability that a given genotype results in a given phenotype. Eq. (<ref>) asserts that the likelihood function p( R | σ, μ) is the Boltzmann distribution of structure R conditional on the amino acid sequence σ. This setting of the likelihood functions p( R | σ, μ) is based on Anfinsen's dogma<cit.>, which states that the state in which a protein adopts its native structure is a thermodynamic equilibrium determined by its amino acid sequence under physiological conditions.

The prior in Eq. (<ref>) is highly non-trivial. To briefly explain the background of our setting of the prior in Eq. (<ref>), it is based on the free energy of a protein with amino acid sequence σ,

    F(σ; β_p, μ_p) = -1/β_plog Z(σ; β_p, μ_p)

which is proportional to the structural partition function Z(σ ; β_p, μ_p) included in it. Since low free energy implies a large partition function, the prior in Eq. (<ref>) can be interpreted as the hypothesis that "amino acid sequences with lower free energy under specific temperature β_p and chemical potential μ_p have evolved." We call this the hypothesis of sequence weights (HSW). HSW was first proposed in <cit.>. Amino acid sequences with high free energy Eq.(<ref>) tend to increase the Hamiltonian Eq. (<ref>) for many structures. The probability that such amino acid sequences form specific compact three-dimensional structures is extremely low. HSW is a hypothesis that preemptively excludes such sequences.

HSW is still an unverified hypothesis, but our previous studies<cit.> have shown that a protein design method assuming HSW exhibits high performance for the 2D lattice HP model with N ≤ 36, which we will analyze in this study. Protein design problem is the inverse problem of protein structure prediction. Additionally, using the prior Eq. (<ref>) based on HSW allows us to cancel out the two partition functions Z(σ ; β, μ) and Ξ(β, μ) required for deriving the posterior, thereby avoiding the computational explosion associated with structural searches.

If β = β_p, μ = μ_p holds, the derivation of the posterior p(σ |  R, μ) is then,

    p(σ| R, μ)    = p( R|σ, μ)p(σ | μ)/∑_σp( R|σ, μ)p(σ|μ)
    
    			   =
    			e^-β H(R , σ ; μ)/Z(σ; β, μ)·Z(σ;β, μ)/Ξ(β, μ)/∑_σe^-β H(R , σ ; μ)/Z(σ;β, μ)·Z(σ;βμ)/Ξ(β, μ)
       = 
                            e^-β H(R , σ ; μ)/Ξ(β, μ)
                            /∑_σ
                            e^-β H(R , σ ; μ)/Ξ(β, μ)
                            
       = e^-β H(R , σ ; μ)/∑_σ e^-β H(R , σ;μ).

From Eq. (<ref>) to Eq. (<ref>), We used that the grand partition function Ξ(β, μ) does not depend on the amino acid sequence σ, allowing it to be factored out of the sum ∑_σ. In the posterior p(σ |  R, μ), the sequence space considered can be regarded, from an evolutionary perspective, as the set of typical sequences that realize a given structure R.

Eq. (<ref>) represents the free energy of a protein with a given amino acid sequence σ. Hence, it is natural to assume that the inverse temperature and chemical potential in this equation are equal to those in the likelihood function Eq. (<ref>), which represents the probability that the protein folds into a given structure R.



 §.§ A free energy depends on a protein structure and its stability to the environmental perturbation



In this subsection, we define the free energy for a sequence space dependent on a particular protein structure R, and derive expressions demonstrating its stability. To facilitate this, we present the expression for the marginal likelihood p( R | μ) here. In Bayesian learning, the marginal likelihood represents the probability of observing the data given a value for the hyperparameter. For the lattice protein case under discussion, the marginal likelihood is

    p( R | μ)    = ∑_σp( R|σ; μ)p(σ | μ) 
       =
    			∑_σe^-β H(R , σ ; μ)/Z(σ; β, μ)·Z(σ; β_p, μ_p)/Ξ(β_p, μ_p)
       =
    			Y( R; β, μ)/∑_ RY( R; β, μ),

where β = β_p and μ = μ_p hold. From Eq. (<ref>), the marginal likelihood, i.e., the probability of observing a protein structure R given the environmental conditions represented by the chemical potential of the surrounding water μ, is proportional to the sequence partition function Y( R; β, μ). Furthermore, the marginal likelihood expressed in Eq. (<ref>) is identical to the denominator of the posterior distribution p(σ| R, μ) when Ξ(β, μ) is not canceled out in the transition from Eq. (<ref>) to Eq. (<ref>). Therefore, the marginal likelihood p( R | μ) serves as the partition function for the posterior distribution p(σ| R, μ) as described in Eq. (<ref>). Consequently, we can consider the corresponding free energy as follows:

    F( R, μ) = -1/βlog p( R | μ).

When a specific structure R is determined at a given inverse temperature β, the free energy F( R, μ) becomes a function solely of the environmental conditions μ. When F( R, μ) is minimized at a particular environmental condition μ = μ_EB—a condition equivalent to maximizing the marginal likelihood, which is referred to as Empirical Bayes estimation in the field of Bayesian learning—the stability of F( R, μ) around μ = μ_EB requires

    ∂^2/∂μ^2 F( R, μ) |_μ=μ_EB > 0.

Before we manipulate Eq. (<ref>), we define the two types of expected values appearing in the transformation of Eq. (<ref>). These are the average taken over the posterior distribution when β = β_p and μ = μ_p hold, and the average taken over the joint distribution, which is the product of the likelihood function and the prior. The joint distribution is given as follows:

    p( R, σ | μ)    = p( R | σ, μ)p(σ | μ) 
       = e^-β H( R, σ;μ)/∑_ R∑_σ e^-β  H( R , σ;μ).

Thus, the joint distribution p( R, σ | μ) forms a Boltzmann distribution where both the structure R and the sequence σ serve as thermal variables. For any physical quantity X( R, σ) that is a function of the structure R and the sequence σ, the average taken over the posterior distribution p(σ |  R, μ) and the average taken over the joint distribution p( R, σ | μ) are given as

    ⟨ X( R, σ) ⟩_|  R   := ∑_σ X( R, σ) p(σ| R, μ) 
       = ∑_σ X( R, σ) e^-β  H( R , σ;μ)/∑_σ e^-β  H( R , σ;μ), 
    ⟨ X( R, σ) ⟩   := ∑_ R∑_σ X( R, σ) p( R, σ| μ) 
       = ∑_ R∑_σ X( R, σ) e^-β  H( R , σ;μ)/∑_ R∑_σ e^-β  H( R , σ;μ).

The notation used on the far left side of Eq. (<ref>) explicitly indicates that the quantity depends on a specific structure R. It is important to note that this does not represent an average over all possible R patterns.

To manipulate inequality (<ref>), we substitute Eq. (<ref>) into Eq. (<ref>), and then transform Eq. (<ref>) as


    < (β∑_i=1^N(1-σ_i))^2>_|  R - (β< ∑_i=1^N( 1 - σ_i)>_|  R)^2   <
          < ( β∑_i=1^N(1-σ_i))^2> - (< β∑_i=1^N( 1 - σ_i)>)^2.


Thus, Eq. (<ref>) results in an inequality between the variance of the number of hydrophilic amino acid residues calculated from the posterior and the variance calculated from the joint distribution. Furthermore, we obtain following expression of Eq. (<ref>):

    β∑_i,j[  < σ_iσ_j >_|  R -  <σ_i >_|  R< σ_j >_|  R]
         <
          β∑_i,j[  < σ_iσ_j > -  <σ_i > < σ_j > ],

where we denote ∑_i,j = ∑_i = 1^N∑_j=1^N. The definition of the magnetic susceptibility χ at the equilibrium state using the spin correlation function is given by χ = β∑_i,j[ <σ_iσ_j> - <σ_i> <σ_j>], where σ_i denotes the i-th spin variable, similar to the current context. Thus, under the posterior Eq. (<ref>), which is the Boltzmann distribution of sequences conditioned on the structure R, so to speak the “hydrophobic susceptibility" can be defined as χ_ R := β∑_i,j[  < σ_iσ_j >_|  R -  <σ_i >_|  R< σ_j >_|  R]. For the joint distribution, which is the Boltzmann distribution over both spaces of R and σ, the hydrophobic susceptibility is denoted as χ := β∑_i,j[  < σ_iσ_j > -  <σ_i > < σ_j > ]. Hence, Eq. (<ref>) becomes the inequality between those two hydrophobic susceptibilities:

    χ_ R < χ.

We here consider the evolutionary meaning of χ_ R. It represents the susceptibility of the hydrophobicity of the structure R to perturbations in μ, within the posterior p(σ |  R, μ). Given that changes in the hydrophobic/hydrophilic composition (hydrophobicity) of amino acid residues can lead to alterations in protein structures <cit.>, χ_ R reflects the structural plasticity of protein R to environmental changes via macroscopic shifts in the gene space that facilitates the formation of R. Therefore, we propose to call χ_ R the evolutionary structural plasticity of a protein structure R. The term evolutionary structural plasticity is used to distinguish it from phenotypic plasticity, which refers to an individual's capacity to adapt to environmental changes.

This quantity, χ_ R, represents the plasticity of a genotype associated with a specific structure R to change for subsequent generations. This concept can be understood by following the simple dynamics of Darwinian evolution. Consider a population of organisms where proteins with the structure R exist in a particular generation. In the subsequent generations, assume that amino acid sequences with slight differences in hydrophobic/hydrophilic composition, which incidentally fold into the same structure R, existed either by chance or due to mutations. Subsequently, perturbations in environmental condition μ occur, leading to the evolution of amino acid sequences that adapt to this new environment, consequently inducing changes in the original structure R. If χ_ R is high, it indicates a higher such genotypic change; conversely, a low χ_ R suggests that changes are less likely. Therefore, χ_ R does not merely reflect the susceptibility of the structure R to change within a single generation but rather its plasticity for change in an evolutionary context.

While there exist proteins whose structures remain unchanged despite variations in hydrophobicity, such proteins are a minority among all proteins. Moreover, it can be assumed that even these exceptional proteins would alter their structures if there were significant differences in hydrophobicity.

χ represents the susceptibility of sequences to change under the joint distribution p( R, σ | μ); thus, χ_ R is akin to an average of overall structural patterns. Therefore, Eq. (<ref>) asserts that for the free energy F( R, μ) of a given structure R to be stable against perturbations in μ, the χ_ R of that structure must be lower than the average χ across all structures.

Finally, for structures that satisfy Eq. (<ref>), we define the following quantity as the steepness of the function around the minimum of the free energy F( R, μ), specifically the second derivative with respect to μ: ∂^2/∂μ^2 F( R,μ) |_μ=μ_EB = χ - χ_ R,

    κ_ R =  χ - χ_ R.

This quantity κ_ R represents the stability against perturbations in μ around the minimum state of the free energy F( R, μ), independent of temperature. 

The evolutionary significance of κ_ R is not immediately apparent. However, since χ_ R represents the evolutionary plasticity in response to environmental changes, κ_ R, as a decreasing function of χ_ R, indicates the robustness of R after several generations under environmental perturbations.

It is important to note that even if hydrophobicity remains constant, different amino acid sequences can lead to different structures. Therefore, κ_ R precisely measures the stability against structural changes due to differences in hydrophobicity, which implies a tolerance for macroscopic structural changes while permitting minor structural variations. The “minor" changes in protein structure that are considered here might include slight alterations due to microscopic differences in the contact network, for example.

For a given structure R, the computation of κ_ R are performed using Belief Propagation (BP) (theoretical details are in Appendix <ref>). BP is an algorithm that efficiently computes the marginal probability of an element in a graph with complex interactions for the entire probability distribution of the graph. The BP algorithm is derived using an extended method of mean-field approximation called the cavity method. It has been shown that the solutions produced by BP are equivalent to the Bethe approximation<cit.>, and if the graph is a tree, the solution is exact. BP also provides good approximate solutions when the graph is close to a tree. Since the contact network of a 2D lattice protein is often a tree graph, BP is suitable for 2D lattice protein models. Indeed, our previous research has shown that the design accuracy of 2D lattice proteins is nearly the same when using BP and Markov Chain Monte Carlo (MCMC) methods<cit.>.



 §.§ Definition of the secondary structure of the
2D lattice proteins



Here we define the secondary structures in 2D lattice proteins, namely α-helices and β-sheets. An α-helix is defined as a structure where amino acid residues form a U-shaped structure connected in alternating orientations and a β-sheet is defined as a structure where amino acid residues are aligned in parallel or anti-parallel configurations (Fig. <ref>). The detection of these secondary structures was performed from contact maps, using the identical method proposed in the previous study <cit.>.


Both α-helices and β-sheets must have at least six residues. This requirement is due to the overrepresentation of single U-shaped structures (four residues) and paired residues in 2D lattice proteins, which would otherwise not qualify as secondary structures. However, a specific case of a 4-residue β-sheet within a 6-residue β-turn (Fig. <ref>c) is included. We count α-helices and β-sheets according to these rules, treating all other regions as random coils.

In lattice proteins, α-helices and β-sheets may share one or two amino acid residues. This study prioritizes β-sheets when two residues are shared and α-helices when one residue is shared. If two residues are shared, the α-helix is resized by removing the shared residues, and if the resized helix has four residues, it is classified as a random coil. If one residue is shared, the β-sheet is resized by removing the shared residue and its paired residue, and if the resized sheet has four residues (unless it is part of a 6-residue β-turn), it is classified as a random coil. Prioritizing β-sheets for two shared residues and α-helices for one shared residue helps balance the two structures.

We do not use 3D lattice proteins because the α-helix in three dimensions requires six residues per turn, unlike the real protein case of 3.6 residues per turn<cit.>. In two dimensions, as shown in Fig. <ref>(a), four residues per turn are closer to realistic proteins. This distinction is crucial as we define the proportion of secondary structures in terms of the number of residues forming these structures within a given protein.

In this study, we use maximally compact structures of N = 6 × 6 residues. The total number of structural patterns is 28,732, excluding structures symmetric under 90^∘ rotation, mirror images along horizontal and vertical axes, diagonal reflections, and head-tail symmetries of the self-avoiding walk.






§ RESULTS



 §.§ Changes in the Number of Lattice Protein Structures According to Evolutionary Stability



First, we examine how the number of lattice protein structures decreases as the stability κ_ R increases, essentially observing how the dimension of the phenotypic space reduces with changes in κ_ R. Specifically, we increment κ_ R from 0 in small steps and record the number of structures that have a κ_ R larger than the value at each step, denoting this number as N^ str(κ_ R). We show the results for the quantity κ_ R divided by β. This measure is taken to cancel out the β in the coefficients of χ_ R and χ, which are common across all structures. We then calculate its proportion out of the total number of structures, which is 28,732. The results, showing the reduction in the number of structures with increasing κ_ R/β within the range 0 ≤κ_ R/β≤ 10, are displayed in Fig.<ref>. In Fig.<ref>, N^tot represents the total number of maximally compact structures for an amino acid residue number N = 6 × 6, equaling 28,732. We set at β =10.


From Fig. <ref>, it is evident that N^ str(κ_ R/β) decreases in a stepwise manner as κ_ R/β increases. It is also observed that there are two points where the rate of decrease in N^ str(κ_ R/β) is particularly significant. Furthermore, approximately 80% of the structures are found within the range 0 ≤κ_ R/β≤ 0.5, indicating that most structures have a gentle slope around the minimum of the free energy F( R, μ). The proportion of structures with κ_ R/β < 0, which are outside the range of Fig. <ref>, is approximately 4% and low. Additionally, there are no structures with κ_ R/β≥ 9.5, as the maximum value of κ_ R/β is 9.3971.




 §.§ Changes in the Proportion of Secondary Structures According to Environmental Robustness



Next, we illustrate the changes in the proportion of secondary structures within groups of lattice protein structures as the stability κ_ R increases incrementally. We also show the results for the quantity κ_ R divided by β. For each increment in κ_ R/β, we plot the average proportions of α-helices, β-sheets, and random coils within the set of structures that have a κ_ R/β greater than the current κ_ R/β value. Fig.<ref> shows its distributions of secondary structures. We also set the inverse temperature at β = 10.


Looking at Fig.<ref>, as κ_ R/β increases, the proportion of α-helices increases while the proportion of β-sheets decreases. Random coils remain neutral to κ_ R/β changes. This result indicates that structures with a higher proportion of α-helices tend to have higher κ_ R/β values and those with a higher proportion of β-sheets tend to have lower κ_ R/β values. Therefore, it can be concluded that proteins with a higher proportion of α-helices are more robust over generations in response to perturbations in environmental conditions.



§ DISCUSSION AND CONCLUSION



First, we note some overall considerations regarding our statistical mechanics theory framed within Bayesian learning for evolution. Our proposed evolutionary theory does not optimize phenotypes to given environmental conditions. That is because the stability analysis around the minimum of the free energy F( R, μ) pertains to the stability around the minimum for environmental conditions μ given a structure (phenotype) R and not the reverse. According to the neutral theory of molecular evolution<cit.>, the evolution of organisms (molecular evolution) is not necessarily a product of optimization under given environmental conditions. Our theory does not contradict these important evolutionary concepts.

The evolutionary significance of κ_ R is, directly, the stability to change in the amino acid sequence σ (genotype) that forms the structure (phenotype) R over subsequent generations. Thus, it cannot necessarily be said to explain the robustness to change of a structure over many generations. Furthermore, κ_ R pertains explicitly to the stability against perturbations in environmental conditions μ, and does not directly inform about the stability against significant changes in μ. Therefore, while κ_ R may be better suited to explaining incremental microevolution, it may not be as applicable to macroevolution. However, it is not uncommon for macroevolution in organisms to be understood as an accumulation of microevolutions; our theory and analysis of κ_ R could still play a significant role in understanding large evolutionary changes.

κ_ R can also be interpreted as part of a broader concept of robustness that includes environmental robustness, especially considering the functionality of gene regulatory networks. In this sense, proposing κ_ R in an evolutionary context is meaningful. However, it should be noted, as emphasized in Chapter II, that all the measures proposed in this study, such as χ_ R and κ_ R, averages over the amino acid sequences σ. Thus, the effects of mutational robustness across successive generations on the structure R remain unclear. This study exclusively analyzed the effects of environmental fluctuations on evolution without integrating comprehensive evolutionary theories that include mutations, natural selection, or genetic drift, which are changes in genes unrelated to natural selection. To include theories that account for mutations and genetic drift, the former would require discussing the stability of quantities that retain amino acid sequence pattern dependencies without integrating over amino acid sequences for proteins. For the latter, it would be necessary to consider quantities or models where hydrophobicity fluctuates solely due to random effects while being constant to changes in the environmental condition μ. Constructing such theories remains a task for future research.

The result of Fig.<ref> indicates that increasing of κ_ R leads to a very rapid reduction in the phenotypic space. It is unclear what proportion of the total structural space, including evolutionary non-viable protein structure patterns, is occupied by actual protein structures, so whether the results from Fig.<ref> accurately explain the dimensional reduction of phenotypes from the perspective of natural selection remains uncertain. However, the fact that evolved phenotypes constitute only a portion of the entire possible phenotypic space can indeed be explained by our proposed metric of robustness to change in phenotypes over subsequent generations in response to perturbations in environmental conditions μ.

Similarly, Fig.<ref> lacks direct biological evidence. However, a conventional study elucidated that α-helices exhibit higher mutational robustness than β-sheets <cit.>. Of course, the relationship between mutational robustness and our results is unclear. However, the fact that α-helices and β-sheets behave differently regarding structural mutational robustness, which is evolutionarily significant, is likely essential. Given the observed variation in the proportions of α-helices and β-sheets among proteins, it is plausible to suggest that these differences may be associated with differences in protein function. Considering that protein function closely relates to evolvability, the differing dependencies on the proportions of α-helices and β-sheets suggest that our proposed stability measure, κ_ R, is likely to be evolutionarily meaningful.

Here, we give the conclusion of our study. We proposed an evolutionary generation model for protein structures by statistical mechanics based on Bayesian learning framework. We considered the chemical potential surrounding proteins as an environmental condition. We discussed the stability of the free energy as a function of protein structure and environmental conditions, as defined in Eq. (<ref>). This stability refers to the robustness to change in hydrophobicity (the proportion of hydrophobic amino acids, analogous to magnetization in magnetic materials), determined by the amino acid sequences that design (fold into) a specific protein structure with high probability. Since changes in hydrophobicity can lead to structural changes, this stability can be considered as the stability of a given protein structure against environmental perturbations over subsequent generations. Consequently, in a 2D square lattice protein model composed of 36 residues, we found that structures with a certain level of stability in their free energy are very rare in the entire structural space, with a higher proportion of α-helices and a lower proportion of β-sheets.



§ ACKNOWLEDGMENTS

This work was supported by JSPS KAKENHI Grant Number 23K19996 (T. T.). The authors are grateful to Macoto Kikuchi, Ayaka Sakata, and Takashi Takahashi for illuminating discussions and helpful comments. This work was also supported by JSPS KAKENHI Grant Number 22H00406 (G. C.), and JSPS KAKENHI Grant Number 22H05117 and JST CREST JPMJCR1912 (Y. K.)




§ THE CALCULATIONS OF HYDROPHOBICITY AND SUSCEPTIBILITY USING BELIEF PROPAGATION









 §.§ The calculations of ⟨σ_i⟩_R and ⟨σ_i⟩


In order to calculate the hydrophobicity of single structure, one needs to obtain the posterior average of single residue <σ_i>_ R given by

    <σ_i>_ R = ∑_σσ_ip(σ| R, μ).

If the posterior p(σ| R, μ) is decoupled to each residue expressed as follows by using the set of residues without σ_i, σ\ i

    p_i(σ_i| R, μ) = ∑_σ\ i p(σ| R, μ),

then, one can rewrite Eq. (<ref>) very simple form as follows:

    <σ_i>_ R   = ∑_σ_i = 0,1σ_ip_i(σ_i| R, μ)
        = p_i(σ_i = 1| R, μ).


Belief propagation (BP) can obtain the marginal distribution p_i(σ_i| R, μ) by using following update rules, 

    ν̃_a → i^(t) (σ_i) = 1/Z_a → i∑_σ_j(i) e^βσ_iσ_j(i)ν_j(i) → a^(t) (σ_j(i)), 
    ν_i → a^(t + 1) (σ_i) = 1/Z_i → a e^βμ(1 - σ_i)∏_b ∈∂_i∖ aν̃_b → i^(t) (σ_i).

In Eqs. (<ref>) and (<ref>), the beliefs or messages ν̃_a → i^(t) (σ_i) and ν_i → a^(t + 1) (σ_i) are the probability from the a-th contact to the i-th residue and the probability from the i-th residue to the a-th contact, respectively. The subscripts a, b, ⋯ are indices on contacts, and the upper right subscript is the number of steps in the BP algorithm. The symbol ∂_i denotes the index set of contacts related to residue σ_i. The constants Z_a → i and Z_i → a are the normalizing constants of each distribution function. The residue index j(i) denotes the index that contacts with i-th residue. In the lattice HP model, all residue-residue interactions are two-body. Thus, this index j(i) is unique to i.

The derivation of the BP update rules (<ref>) and (<ref>) are somewhat technical, so they are not presented here. Please refer to the appendix in our previous work<cit.> for information of the derivation of the above BP update rules.

If one properly defines ν_i → a^(t = 0)(σ_i) as the initial condition (in this study, we use 0.5.) and computes Eqs. (<ref>) and (<ref>) at each step for all combinations (i, a), after sufficient iterations t_ max, the following belief:

    ν_i(σ_i) = 1/Z_i∏_a ∈∂_iν̃_a → i^(t_ max) (σ_i),

converges to the marginal distribution p_i(σ_i |  R,μ). In Eq. (<ref>) where Z_i is the normalization constant.

In order to obtain another hydrophobicity h, one needs to calculate the joint average of σ_i. It is expressed as follows.

    <σ_i>    = ∑_ R∑_σσ_ip( R, σ | μ)
        = ∑_ R∑_σσ_ip( R | μ)p(σ |  R, μ)
        = ∑_ R∑_σσ_ip(σ |  R, μ)p( R | μ)
        = ∑_ R<σ_i>_ R p( R | μ)
        = ∑_ R<σ_i>_ R Y( R;β,μ)/∑_ RY( R;β,μ).

In Eq. (<ref>) to Eq. (<ref>), we used Eq. (<ref>). Therefore, one has to calculate the sequence partition function Y( R;β,μ) = ∑_σe^-β H( R, σ; μ) and the posterior average of each residue <σ_i>_ R for all lattice structure patterns. We have the all patterns of the self avoiding walks on 2D N = 6×6 square. Thus, in this study, we can obtain the exact value of < σ_i>.

For more large size of lattice proteins, one has to carry out the efficient multi-canonical Monte Carlo methods suitable for exploring lattice protein structures<cit.>. For the realistic protein structures, such a structural search is extremely difficult even with the use of the database<cit.>, because the structural search space has to involve the random structural patterns including structures that have not evolved.

By using BP, The sequence partition function Y( R;β,μ) = ∑_σe^-β H( R, σ; μ) is obtained from the Bethe free entropy F_B(ν̃^*) = log Y( R;β,μ) where ν̃^* is the set of the contact to residue messages (Eq. (<ref>)) after sufficient time steps. Free entropy is -β times free energy.

The Bethe free entropy F_B(ν̃^*) is the free entropy under the Bethe-approximation. In the lattice HP model using the Hamiltonian (1), F_B(ν̃^*) is given by

    F_B(ν̃^*)    = ∑_a=1^Mlog Z_a + ∑_i=1^Nlog Z_i - ∑_ialog Z_ia,
    
        Z_a   := ∑_σ_iσ_j(i) e^βσ_iσ_j(i)∏_i ∈∂ aν_i → a^*(σ_i),
    
        Z_i   := ∑_σ_i e^βμ (1-σ_i)∏_b ∈∂ iν̃_b → i^*(σ_i),
    
        Z_ia   := ∑_σ_iν̃_a → i^*(σ_i) ν_i → a^*(σ_i),

where ia denotes index of elements of the set in which all combinations of residues and contacts, and the symbol ∂ a denotes the index set of residues related to contact a. The symbol M denotes the number of contacts. The messages ν_i → a^*(σ_i) and ν_a → i^*(σ_i) are the converged i-th residue to a-th contact message and the converged a-th contact to i-th residue message, respectively. The derivation of Eqs. (<ref>) ~(<ref>) for general cases involves a technical and lengthy explanation, which is beyond the scope of this paper. For further details, please refer to the representative texts<cit.>. Then, we obtain

    Y( R;β,μ) = ( ∏_a Z_a) ( ∏_i Z_i)/∏_ia Z_ia.


In order to obtain μ_EB for each structure, we derive the minimization condition of free energy of a structure F( R, μ) defined by Eq. (<ref>). The minimization condition with respect to μ : (∂ / ∂μ) F( R, μ) = 0 becomes as follows:

    < ∑_i=1^Nσ_i>_ R = < ∑_i=1^Nσ_i>.

The minimization parameter μ_EB satisfies Eq. (<ref>). Thus, one can get μ_EB by computing ∑_i< σ_i>_ R and ∑_i< σ_i> through the method explained above for each structure. We used the bisection method to compute left hand side and right hand side of Eq. (<ref>).



 §.§ The calculations of χ_R and χ


In order to obtain the susceptibility χ_ R, one has to calculate the residue correlation function < σ_iσ_j>_ R given by

    <σ_iσ_j>_ R = ∑_σσ_iσ_jp(σ| R, μ).

If σ_i and σ_j are not statistically independent, i.e., when there is at least one path connecting σ_i and σ_j on the contact graph, one has to think the joint probability

    p_i,j(σ_i, σ_j |  R, μ) = p_i | j(σ_i | σ_j,  R, μ) p_j(σ_j |  R, μ),

where the conditional probability, p_i | j(σ_i | σ_j,  R, μ) is given by 

    p_i | j(σ_i | σ_j,  R, μ) = p_i,j(σ_i, σ_j |  R, μ)/p_j(σ_j |  R, μ).

If the following marginalization can be carried out

    p_i | j(σ_i | σ_j,  R, μ) = ∑_σ\ ip(σ | σ_j , R,μ),

one can then obtain < σ_iσ_i>_ R following quite simple form

    < σ_iσ_i>_ R   = ∑_σ_i∑_σ_jσ_iσ_j p_i | j(σ_i | σ_j,  R, μ) p_j(σ_j |  R, μ)
       = p_i | j(σ_i = 1 | σ_j=1,  R, μ) p_j(σ_j=1 |  R, μ).

BP is also able to calculate the conditional marginal p_i | j(σ_i | σ_j,  R, μ). The procedure is as follows: one uses the conditional messages ν̃_a → i|j^(t) (σ_i | σ_j = σ) and ν_i|j → a^(t + 1) (σ_i | σ_j = σ) instead of normal BP messages Eqs.(<ref>) and (<ref>), where σ is the realization of σ_j. Thus, in the current case σ = 1, the normal BP messages change to following conditional forms

    ν̃_a → i|j^(t) (σ_i | σ_j = 1)   = 1/Z_a → i|j∑_σ_k(i) e^βσ_iσ_k(i)ν_k(i)|j → a^(t) (σ_k(i)|σ_j = 1),
    ν_i|j → a^(t + 1) (σ_i | σ_j = 1)   = 1/Z_i|j → a e^βμ(1 - σ_i)∏_b ∈∂_i∖ aν̃_b |j→ i^(t) (σ_i | σ_j = 1).

In Eq. (<ref>), we use k(i) as the index of residue that contact with σ_i to avoid confusion with j of current context. The symbol Z_a → i|j and Z_i|j → a are the  normalization constants for the corresponding messages. As with the normal BP described earlier, the following belief:

    ν_i|j(σ_i | σ_j = 1) = 1/Z_i|j∏_a ∈∂_iν̃_a → i|j^(t_ max) (σ_i | σ_j = 1),

converges to the conditional marginal p_i | j(σ_i | σ_j=1,  R, μ). The symbol Z_i|j in Eq. (<ref>) denotes the normalization constant of this message.

One can calculate χ by using following formula

    < σ_iσ_j> = ∑_ R<σ_iσ_j>_ R Y( R;β,μ)/∑_ RY( R;β,μ).

Eq. (<ref>) is obtained by the same manner as the derivation process of Eq.(<ref>).



















































































apsrev4-2



	





