Disentanglement of Evolutionary Constraints in Statistical Models of Proteins

The exponential growth of protein sequences in the post-genomic era has revolutionized the application of generative sequence models for pivotal tasks such as contact prediction, protein design, alignment, and homology search. Despite remarkable progress in these areas, the interpretability of the modeled pairwise parameters remains limited due to complexities arising from coevolution, phylogeny, and entropy. While post-correction methods for contact prediction have been developed to eliminate entropy-related contributions from predicted contact maps, there is currently no direct approach to correct entropy in other applications reliant on raw parameters. In this paper, we investigate the sources of entropy signal and propose a novel spectral regularizer, LH (an abbreviation of Henri Lebesgue), to mitigate its impact during model fitting. By incorporating this regularizer into the GREMLIN framework (utilizing a Markov random field or Potts model), we enable the accurate inference of sparse contact maps while simultaneously improving interpretability and addressing overfitting concerns critical for sequence evaluation and design. To validate the efficacy of our approach, we design multiple protein sequences based on GREMLIN with both L2 and LH regularizers, and subsequently experimentally measure their using cDNA display proteolysis. Our findings demonstrate that proteins designed using the LH regularizer exhibit increased diversity and enhanced folding stability.

Figure 1

In addressing the entropy issue, it is essential to clarify that here, entropy pertains to the single-site entropy within each column of the multiple sequence alignment (MSA). The MSA is utilized as input for the Markov random field (MRF) model, aiming to reconstruct the original sequences. The model comprises a tensor representing coevolution ($L\times K\times L\times K$) and a matrix capturing conservation and positional entropy ($L\times K$), where $L$ is the length of the protein sequence, and $K$ is the number of amino acids plus an aligned gap. Typically, the four-dimensional tensor is condensed into an $L\times L$ matrix, followed by the application of average product correction (APC) to disentangle contacts from the entropy signal. Traditionally, in applications like the critical assessment of structure prediction (CASP13), a prominent blind protein structure prediction experiment, the reduced $L\times L$ matrix is favored over the full four-dimensional (4D) tensor as input for deep learning. This coevolution model serves not only for structural predictions but also for predicting mutation effects and aligning sequences. However, it is crucial to note that the entropy signal (highlighted in red) is exclusively separated through the APC method at the reduced matrix level. At the 4D tensor level, the coevolution signal and entropy remain intertwined, posing potential challenges to various applications, as illustrated by the dashed red lines. This entanglement at the 4D tensor level may impact the overall performance of related applications.

Figure 2

The effects of regularization type (L2, LH, and LB) on the sparsity of contact maps, gradient, and overfitting. The first row [(a)–(c)] shows the $\mathbf{M}$, with L2, LH, and LB regularization, respectively. The second row [(d)–(f)] shows the gradient of these three regularizers. Gradients are consistently recorded throughout the entirety of the training process, and the depicted gradients correspond to those obtained from the final training step as it approaches convergence. The third row [(g)–(i)] shows the distributions of reconstruction losses for the natural MSA training set (blue), the test set (yellow), MSA sampled from a PSSM of the natural MSA (green), and MSA sampled from a random distribution MSA (red). Even though the loss distribution of the training set is well separated from the PSSM and randomly sampled sequences, for the L2 and LB regularized models the test set distribution does not overlap with the training set loss, indicating an overfitting issue. For the LH regularized model, the training and test set loss distributions have a good overlap.

Figure 3

LH regularized parameters disentangle conservation and entropy from coevolution while L2 regularized parameters do not. The figure shows the distribution of Pearson correlation values across different proteins, comparing the PSSMs of the sampled sequences vs the natural sequences. The sampling was done using just one-body parameters (b), just two-body parameters (w), or a combination of both one- and two-body parameters (wb).

Figure 4

For unsupervised protein contact prediction, the average product correction (APC) is no longer required under the LH regularizer. (a) The accuracy of the L2-based method is set as a baseline with two solid lines. The blue one is the performance without the APC, and the orange one is the accuracy with the APC. The $x$-axis is the hyperparameter for LH, and the red dots showed the raw accuracy of LH and blue dots with APC correction. (b) The accuracy of Block L1 is shown with different hyperparameter scanning. (c) The best raw performance from LH compared with the L2-based method. Each point corresponds to a protein MSA, the axes indicate the accuracy of each method, defined as the average precision of the top $L$-ranked contacts, and $L$ is the length of the protein. (d) Performance comparison between L2 and LB.

Figure 5

For protein contact prediction, three different ways of contact extraction are applied to all three methods, L2 (green), LH (red), and LB (yellow). Besides the raw and APC method, we treated the $\mathbf{W}$ matrices as inputs of a simple regression model. The top $L$ precision is used to curve the performance, and $L$ is the length of the protein.

Figure 6

Reanalysis of Spearman correlation between the Hamiltonian of GA/GB and published folding temperature [19]. The L2 and LH regularizers are applied to curve the Hamiltonian. The $x$-axis is the Hamiltonian, and the $y$-axis is the reported Tm. Each point is a designed protein sequence.

Figure 7

The histogram of $\mathrm{\Delta }G$ difference between LH and L2 with different optimization protocols. The white line represents a scaling adjustment for the illustration. The red line corresponds to the natural protein folding rate and serves as a designed baseline for comparison. It shows that, no matter which optimization protocol is adopted, $\mathrm{\Delta }G$ of proteins designed by LH are larger than those designed by L2, in most of the cases.

Figure 8

The association between $\mathrm{\Delta }\mathrm{\Delta }G$ and sequence identity is examined across proteins designed using various models and optimization protocols. The analysis reveals that proteins designed with the LH model exhibit greater sequence diversity and stability compared to those designed with the L2 model. Instances where $\mathrm{\Delta }\mathrm{\Delta }G$ values surpass the dashed lines signify that the designed proteins are more stable than naturally occurring proteins.

Figure 9

PCA analysis of MSAs to check about the sequence relations. All sequences are projected into PC1 and PC2; the first plot is colored by norm r.e. (fitness), and the rest are colored by a negative Hamiltonian. Blue couples with positive fitness.

Figure 10

Structural reanalysis [42] of thermal stability from (a) ANC1 to ANC4, (b) ANC1 to B. marinas, and (c) ANC1 to B. subtilis using published crystal structures [42, 44, 45]. The right panel showed the difference in the coevolutionary coupling strength of two different proteins. The more unstable the right pairwise interaction is, the more blue the pattern appears. The solid red lines represent high thermal stability, and the dashed-dotted line weakened the interactions for the right protein.