# Generative Capacity of Probabilistic Protein Sequence Models

Francisco McGee  Center for Biophysics and Computational Biology, Temple
University, Philadelphia, 19122, USA  Institute for Computational Molecular
Science, Temple University, Philadelphia, 19122, USA  Department of Chemistry,
Temple University, Philadelphia, 19122, USA  Quentin Novinger  Institute for
Computational Molecular Science, Temple University, Philadelphia, 19122, USA
Department of Computer & Information Sciences, Temple University,
Philadelphia, 19122, USA  Ronald M Levy  Center for Biophysics and
Computational Biology, Temple University, Philadelphia, 19122, USA  Department
of Biology, Temple University, Philadelphia, 19122, USA  Department of
Chemistry, Temple University, Philadelphia, 19122, USA  Department of Physics,
Temple University, Philadelphia, 19122, USA  Vincenzo Carnevale  Institute for
Computational Molecular Science, Temple University, Philadelphia, 19122, USA
Department of Biology, Temple University, Philadelphia, 19122, USA
Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS>Allan Haldane  Center for Biophysics and Computational Biology, Temple
University, Philadelphia, 19122, USA  Department of Physics, Temple
University, Philadelphia, 19122, USA  Corresponding authors:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Potts models and variational autoencoders (VAEs) have recently gained
popularity as generative protein sequence models (GPSMs) to explore fitness
landscapes and predict the effect of mutations. Despite encouraging results,
quantitative characterization and comparison of GPSM-generated probability
distributions is still lacking. It is currently unclear whether GPSMs can
faithfully reproduce the complex multi-residue mutation patterns observed in
natural sequences arising due to epistasis. We develop a set of sequence
statistics to assess the “generative capacity” of three GPSMs of recent
interest: the pairwise Potts Hamiltonian, the VAE, and the site-independent
model, using natural and synthetic datasets. We show that the generative
capacity of the Potts Hamiltonian model is the largest; the higher order
mutational statistics generated by the model agree with those observed for
natural sequences. In contrast, we show that the VAE’s generative capacity
lies between the pairwise Potts and site-independent models. Importantly, our
work measures GPSM generative capacity in terms of higher-order sequence
covariation statistics which we have developed, and provides a new framework
for evaluating and interpreting GPSM accuracy that emphasizes the role of
epistasis.

## Introduction

Recent progress in decoding the patterns of mutations in protein multiple
sequence alignments (MSAs) has highlighted the importance of mutational
covariation in determining protein function, conformation and evolution, and
has found practical applications in protein design, drug design, drug
resistance prediction, and classification [1, 2, 3]. These developments were
sparked by the recognition that the pairwise covariation of mutations observed
in large MSAs of evolutionarily diverged sequences belonging to a common
protein family can be used to fit maximum entropy “Potts” statistical models
[4, 5, 6]. These contain pairwise statistical interaction parameters
reflecting epistasis [7] between pairs of positions. Such models have been
shown to accurately predict physical contacts in protein structure [6, 8, 9,
10], and have been used to significantly improve the prediction of the fitness
effect of mutations to a sequence compared to site-independent sequence
variation models which do not account for covariation [11, 12]. They are
“generative” in the sense that they define the probability, $p(S)$, that a
protein sequence $S$ results from the evolutionary process. Intriguingly, the
probability distribution $p(S)$ can be used to sample unobserved, and yet
viable, artificial sequences. In practice, the model distribution $p(S)$
depends on parameters that are found by maximizing a suitably defined
likelihood function on observations provided by the MSA of a target protein
family. As long as the model is well specified and generalizes from the
training MSA, it can then be used to generate new sequences, and thus a new
MSA whose statistics should match those of the original target protein family.
We refer to probabilistic models that create new protein sequences in this way
as generative protein sequence models (GPSMs).

The fact that Potts maximum entropy models are limited to pairwise epistatic
interaction terms and have a simple functional form for $p(S)$ raises the
possibility that their functional form is not flexible enough to describe the
data, i.e. that the model is not well specified. While a model with only
pairwise interaction terms can predict complex patterns of covariation
involving three or more positions through chains of pairwise interactions, it
cannot model certain triplet and higher patterns of covariation that require a
model with more than pairwise interaction terms [13]. For example, a Potts
model cannot predict patterns described by an XOR or boolean parity function
in which the $n$-th residue is determined by whether an odd number of the
$n-1$ previous residues have a certain value (see Supplementary Information).
While some evidence has suggested that in the case of protein sequence data
the pairwise model is sufficient and necessary to model sequence variation
[14, 15, 16], some of this evidence is based on averaged properties, and there
appears to be some weak evidence for the possibility of rare “higher-order
epistasis” affecting protein evolution [17, 7, 18, 19], by which we mean the
possibility that subsequence frequencies of three or more positions cannot be
reproduced by a model with only pairwise interactions. Fitting maximum entropy
models with all triplet interactions is not feasible without significant
algorithmic innovation, since for a protein of length 100 it would require
approximately 10B parameters and enormous MSA datasets to overcome finite
sampling error (see Supplementary Information). However, recent developments
in powerful machine learning techniques applied to images, language, and other
data have shown how complex distributions $p(S)$ can be fit with models using
more manageable parameter set sizes. Building on the demonstrated power of
incorporating pairwise epistasis into protein sequence models, this has
motivated investigation of machine learning strategies for generative
modelling of protein sequence variation which can go beyond pairwise
interactions, including Restricted Boltzmann Machines (RBMs) [3], variational
autoencoders (VAEs) [20, 21, 22, 23], General Adversarial Networks (GANs)
[24], transformers [25, 26, 27, 28], and others.

One model in particular, the VAE [29, 30], has been cited as being well suited
for modelling protein sequence covariation, with the potential to detect
higher order epistasis [20, 21]. The VAE also potentially gives additional
insights into the topology of protein sequence space through examination of
the “latent” (hidden) parameters of the model, which have been suggested to be
related to protein sequence phylogenetic relationships [20, 21, 22, 27]. One
implementation of a VAE-GPSM, “DeepSequence”, found that the VAE model was
better able to predict experimental measurements of the effect of mutations
reported in deep mutational scans than a pairwise Potts model, which was
attributed to the VAE’s ability to model higher-order epistasis [20, 11].
However, it has also been suggested by others that the improvement reported
for DeepSequence could be attributed to the use of biologically motivated
priors and engineering efforts, rather than because it truly captured higher-
order epistasis [21]. Furthermore, while VAE-GPSMs are generative and aim to
capture the protein sequence distribution $p(S)$, to our knowledge none of
these studies have thoroughly tested what we will call the “generative
capacity” [31, 32] of the VAE model, meaning the ability of the model to
generate new sequences drawn from the model distribution $p(S)$, which are
statistically indistinguishable from those of a given “target” protein family.
Testing the generative capacity, specifically higher-order covariation, of a
GPSM is a fundamental check of whether the model is well specified and
generalizes from the training set, two prerequisites to capturing higher-order
epistasis.

Here, we perform a series of numerical experiments that explore the generative
capacity of GPSMs. These GPSMs include a pairwise Potts Hamiltonian model with
pairwise interaction terms (Mi3) [33], two state-of-the-art implementations of
a variational autoencoder (VAE), and a site-independent model which does not
model covariation (Indep). One VAE implementation is a “standard" VAE
implementation (sVAE), which uses an architecture very similar to that of
EVOVAE[22] recently used to model protein sequences, though optimized with a
larger number of latent dimensions, and which more closely follows the VAE
inference method as it was originally presented [29, 30]. The second VAE
implementation is “DeepSequence”, mentioned above, which uses a deeper neural
network and a more sophisticated optimizer (see Methods) [20].

We evaluate the generative capacity of these models using three standard MSA
statistics, and introduce a fourth and novel metric based on averaged higher-
order marginals ($r_{20}$) [16]. The three standard metrics are pairwise
covariance correlations [2, 16, 34, 35], Hamming distance distributions [2,
16, 36, 37], and statistical energy correlations [16, 20, 38]. We also test
how each GPSM behaves as fewer training sequences are provided, and thus how
well the GPSMs generalize, using two training MSA sizes: one representing the
estimated size upper bound for a Pfam protein family (10K) (see Supplementary
Information) [39], and one being large enough to eliminate out-of-sample error
in our generative capacity measurements (see Results) [40]. We evaluate GPSM
performance both for natural datasets and for synthetic datasets in which the
ground truth is known. Our comparative analysis thus consists of four
intersecting variables: three GPSMs, a suite of four generative capacity
metrics, two training MSA sizes, and two dataset types.

Our results show that while VAEs model some epistasis they are unable to
reproduce MSA statistics as well as Mi3, but are more accurate than Indep,
which cannot model covariation, for all metrics tested. The two VAE
implementations we tested perform similarly to each other across all metrics.
We find that the VAEs perform more similarly to Mi3 than they do to Indep, but
perform more like the Indep model for the $r_{20}$ metric than does Mi3, which
is the metric that most directly tests a GPSM’s ability to model higher-order
sequence covariation. This result suggests that $r_{20}$ is a more sensitive
measure of GPSM generative capacity than the other standard metrics,
representing a novel and powerful metric to discriminate between GPSMs by
averaging over higher-order covariation terms up to the 10th order. Our work
consolidates several benchmark statistics of generative capacity for GPSMs and
introduces a new one, $r_{20}$, offering a novel framework for evaluating and
interpreting GPSM accuracy in the context of higher-order covariation. By
quantifying and comparing GPSM performance in our innovative epistasis-
oriented approach, we hope to better understand the challenges and limitations
inherent to generative modelling of natural protein sequence datasets, better
gauge the state of the art, and provide insight for future efforts in terms of
minimizing the confounding effects of data limitations in generative protein
sequence modeling.

## Results

### Target Distributions

Our goal is to set baseline expectations for the generative capacity of the
GPSMs when fit to synthetic or natural protein sequence data of varying
training MSA sizes. Generative models of protein MSAs define a distribution
$p_{\theta}(S)$ for the probability of a sequence $S$ appearing in an MSA
dataset given model parameters $\theta$. The model parameters are fit by
either exact or approximate maximum likelihood inference of the likelihood
$\mathcal{L}=\sum_{S\in\text{MSA}}p_{\theta}(S)$ over a training MSA, using
regularization techniques to prevent overfitting. The sequences in the
training MSA are assumed to be identical independent samples from a “target”
probability distribution $p^{0}(S)$, which is generally unknown [29]. For a
model with high generative capacity, $p_{\theta}(S)$ will closely approximate
$p^{0}(S)$ [21]. We test three GPSMs: a pairwise Potts model (Mi3) [33]; two
implementations of a VAE, which are a standard VAE (sVAE) and DeepSequence
[20]; and a site-independent model (Indep), each with a different functional
form of $p_{\theta}(S)$.

It is not possible to measure the similarity of $p_{\theta}(S)$ and $p^{0}(S)$
directly because of the high dimensionality of sequence space, since the
number of sequence probabilities to compare is equal to $q^{L}$, where $L$ is
the sequence length (typically $\sim$300) and $q$ is the alphabet size
($\sim$21). Instead, we measure how derived statistics computed from
evaluation MSAs generated by each GPSM match those of target MSAs drawn from
the target distribution. Three of these are standard metrics: the pairwise
Hamming distance distribution, the pairwise covariance scores, and the GPSM’s
ability to predict $p^{0}(S)$ for individual sequences, also called
statistical energy (see Methods). We present a fourth metric, the averaged
higher order marginal accuracy $r_{20}$.

### GPSM Error and Experiments

Figure 1: 5-stage analysis pipeline diagram. Stage 1: Data Sequestration. Two different MSAs are sequestered: one natural from Uniprot/TREMBL; and the other synthetic, generated by a Potts model fit to the natural MSA. Stage 2: Data Processing. Sequences are indexed and split. Non-overlapping training, target, and evaluation MSAs are shown. Stage 3: GPSM Training. GPSMs are trained. Stage 4: Artificial Sequence Generation. Evaluation MSAs are generated from each GPSM. Stage 5: Generative Capacity Measurements. Computation and visualization of generative capacity metrics is performed. Specification Error |  occurs when the functional form $p_{\theta}(S)$ of a GPSM is not flexible enough to accurately model the target distribution $p^{0}(S)$ for any choice of parameters $\theta$
---|---
Out-of-Sample Error |  occurs when a GPSM fit to a finite training dataset fails to correctly model unseen data, and is a consequence of overfitting
Estimation Error |  occurs due to statistical error in the MSA evaluation metrics when computed from finite evaluation and target MSAs
Training MSA |  drawn from $p^{0}(S)$, used to train and parameterize a GPSM
Target MSA |  drawn from $p^{0}(S)$ separately from the training MSA, used as a validation dataset for performing metrics
Evaluation MSA |  drawn from $p_{\theta}(S)$ for a parameterized GPSM, used to compare to the target MSA
Table 1: Glossary of error types, and the MSA datasets we use to evaluate
these errors.

We divide our tests into two analysis tracks: one synthetic, in which we train
the GPSMs on synthetic MSAs generated from a known target distribution; and
one natural, in which we train the GPSMs on a representative natural protein
family MSA, the kinase super family, sequestered from Uniprot/TrEMBL [41].
These analysis tracks are meant to probe and isolate two distinct forms of
error which may cause $p_{\theta}(S)$ to deviate from $p^{0}(S)$ (see Table
1). The first is “specification error” [42], which occurs when the functional
form of $p_{\theta}(S)$ of a model is not flexible enough to accurately model
the target probability distribution $p^{0}(S)$ for any choice of parameters. A
key motivation for choosing a VAE over a Potts model is its potentially lower
specification error when higher-order epistasis is present [20]. Indeed, Potts
models are limited to pairwise interaction terms of a particular functional
form, while VAEs are not. The second form of error is “out-of-sample error”
[43], caused by a paucity of training samples, and is the consequence of
overfitting [44]. Even a well-specified model could fail to generalize when
fit to a small training dataset, and may mis-predict $p^{0}(S)$ for test
sequences, so it follows that increasing training MSA size reduces out-of-
sample error. Beyond specification and out-of-sample error, which each reflect
an aspect of GPSM generalization error, there can be “estimation error” in our
MSA test statistics due to the finite MSA sizes we use to estimate their
values, which sets an upper bound on how well these statistics can match their
target values, depending on the metric [45]. Finally, other errors may arise
due to implementation limitations of the inference methods, for instance due
to finite precision arithmetic or to finite sample effects when Monte Carlo
methods are used.

The synthetic analysis allows us to isolate specification error, and minimize
both out-of-sample and estimation error, because here the target distribution
is known exactly, and we can generate arbitrarily large training, target, and
evaluation datasets. The synthetic analysis also allows us to quantify out-of-
sample error by modulating the training MSA size. We specify the synthetic
target probability distribution $p^{0}(S)$ to be exactly the Potts model
distribution we inferred based on natural protein kinase sequence data using
Mi3 in our natural analysis (see Methods) [33]. The sequences generated from
this synthetic target distribution should have statistical properties similar
to real, or “natural”, protein family MSAs, albeit constrained by the fact
that the Hamiltonian model used to generate the synthetic dataset is limited
to pairwise epistatic interaction terms only [6, 16, 46]. Whereas Indep and
the VAEs may still fail to model this known probability distribution in the
synthetic analysis, our expectation is that Mi3 will be unaffected by
specification error in the synthetic tests, since the target MSA is sampled
from the same probability distribution used to carry out the inference. In
Supplementary Information, we also perform an alternate synthetic test which
does not favor Mi3 in this way, in which the target distribution is instead
specified by sVAE, finding that both Mi3 and sVAE are able to fit this target
sVAE distribution accurately. In our synthetic analysis, we test two synthetic
training MSA sizes: (i) 1M sequences, to minimize overfitting effects and
consequently out-of-sample error, thereby isolating GPSM specification error
in this experiment; and (ii) 10K sequences, to illustrate the expected GPSM
performance on typical datasets, as most protein families in Pfam have less
than 10K independent effective sequences (see Supplementary Information) [39].

The natural analysis examines the performance of the models on natural
sequence data, which potentially contain higher-order correlations that
require a Hamiltonian model with triplet or higher-order interaction terms to
capture. On this dataset, the VAEs could potentially outperform Mi3, depending
on the importance of higher-order epistatic terms, if present [20]. However,
unlike in the synthetic analysis, here we do not know a priori the target
distribution and, most importantly, we have only limited datasets for both
training and evaluation. In the natural tests, the training and target MSAs
each contain $\sim$10K non-overlapping kinase sequences from the
Uniprot/TREMBL database after phylogenetic filtering at 50% sequence identity.
After this filtering we consider each sequence to be an independent sample of
the evolutionary process, which has the equilibrium distribution $p^{0}(S)$.

Our overall testing procedure is outlined in Fig.1 (see Methods), and the
terms discussed are summarized in Table 1. Our training datasets are either a
natural protein sequence dataset obtained from Uniprot/TREMBL, or a synthetic
training dataset. We then fit the GPSMs to the training datasets, and generate
evaluation MSAs from each model. Finally, using our suite of four generative
capacity metrics, we compare statistics of the evaluation MSAs to those of
“target” MSAs, which contain sequences drawn from the target distributions and
therefore represent our expectation.

### Pairwise covariance correlations

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Pairwise covariance correlations (Top Row) and Pearson $r_{20}$
correlations (Bottom Row). GPSM MSA statistics are compared to the
corresponding target values. GPSMs were trained on 1M synthetic (a, d), 10K
synthetic (b, e), or 10K natural (c, f) kinase sequences. Top Row Covariances
$C^{ij}_{\alpha\beta}$ computed from MSAs of 500K sequences generated by each
GPSM (y-axis) vs target covariances (x-axis) computed from MSAs of 500K
sythnetic target sequences (a, b) or from 10K natural target sequences (c),
with Pearson correlation $\rho$ shown for each comparison. Generative capacity
slightly decreases as synthetic training sample size is reduced from 1M in a
to 10K in b for all GPSMs, except Indep. Covariances around zero are omitted
for simplicity. Bottom Row Pearson $r_{20}$ score (y-axis) as a function of
Higher-Order-Marginal (HOM) $n$-tuple length (x-axis), illustrating how well
each GPSM predicts HOMs from the target distribution. Lines are colored as in
the top panels. $r_{20}$ scores for each GPSM are computed using GPSM-
generated evaluation MSAs of 6M sequences compared to target MSAs of 6M
sequences for the synthetic tests (d, e) or to a target MSA of 10K natural
sequences for the natural test (f) due to limited natural sequence data. The
black dotted line denotes an estimation upper limit caused by finite sampling
estimation error. This is computed as the $r_{20}$ scores for two non-
overlapping synthetic target MSAs of 6M sequences each (d and e, black
circles) or for MSAs of 6M and 10K sequences both generated from the Mi3 model
trained on natural sequences (f, black triangles). For panel f, only lengths
two through seven are plotted, as the small dataset size limits HOM
estimation. Insets emphasize pairwise $r_{20}$. Comparing panels d and e, the
generative capacity of all GPSMs decreases with synthetic training sample size
according to $r_{20}$.

We first examine the pairwise covariance scores for pairs of amino acids of an
MSA defined as
$C^{ij}_{\alpha\beta}=f^{ij}_{\alpha\beta}-f^{i}_{\alpha}f^{j}_{\beta}$. Here,
$f^{ij}_{\alpha\beta}$ are the MSA bivariate marginals, meaning the frequency
of amino acid combination $\alpha,\beta$ at positions $i,j$ in the MSA.
$f^{i}_{\alpha}$ and $f^{j}_{\beta}$ are the univariate marginals, or
individual amino acid frequencies at positions $i$ and $j$. Each covariance
term measures the difference between the joint frequency for pairs of amino
acids and the product of the single-site residue frequencies, i.e. the
expected counts in the hypothesis of statistical independence. The scores
equal 0 if the two positions do not covary. Coevolving amino acids are an
important aspect of sequence variation in protein MSAs, and a GPSM’s ability
to reproduce the pairwise covariance scores of the training dataset has been
used in the past as a fundamental, non-trivial measure of the GPSM’s ability
to model protein sequence covariation [2, 16, 34, 35].

For each GPSM, we compare pairwise covariance scores for all pairs of
positions and residues $\hat{C}^{ij}_{\alpha\beta}$ in their respective
evaluation MSA to the corresponding target pair $C^{ij}_{\alpha\beta}$ in the
target MSA using the Pearson correlation coefficient
$\rho(\\{C^{ij}_{\alpha\beta}\\},\\{\hat{C}^{ij}_{\alpha\beta}\\})$ (Fig. 2,
Top Row). In the synthetic tests we evaluate this statistic using 500K
sequences for both the target and evaluation MSAs, while for the natural test
we compare 500K evaluation sequences to the available 10K target sequences.
Mi3 accurately reproduces the target covariance scores in all tests
($\rho=0.99$ in Fig. 2(a), $\rho=0.95$ in Fig. 2(b), and $\rho=0.88$ in Fig.
2(c) respectively). The somewhat lower value for Mi3 in the natural analysis
of $\rho=0.88$ is accounted for entirely by increased estimation error in that
test, as only 10K target sequences are available for evaluation, and the
expected $\rho$ due only to estimation error is $\rho\sim$0.87, which is
computed by comparing the natural 10K target sequences to the natural 10K
training sequences using this metric. The high generative capacity of Mi3 is
expected, because Mi3 parameters are optimized to exactly reproduce the joint
frequencies of pairs of amino acids from the target MSA. The VAEs’ inference
does not include this constraint, but we find that even when trained on the
larger (1M) dataset of synthetic sequences, the covariances computed from the
VAEs’ evaluation MSAs are generally similar to each other, and smaller in
magnitude than those of the target, showing smaller correlation with the
target than Mi3 ($\rho=0.9$ for sVAE and $\rho=0.92$ for DeepSequence in Fig.
2(a)). For the VAEs, this amount of error in $\rho$ can primarily be
attributed to specification error, since training GPSMs on 1M sequences
largely eliminates out-of-sample error, and the large evaluation MSAs make
estimation error negligible. The VAEs’ covariances are further scaled down
slightly when fit to the synthetic 10K dataset ($\rho=0.87$ for sVAE and
$\rho=0.89$ for DeepSequence in Fig. 2(b)). Indep cannot reproduce covariances
by definition, so $\rho$ is zero in all tests, as expected. The generative
capacity trends for this metric are consistent between the synthetic and
natural analyses for all GPSMs, showing the behavior is not due to artificial
properties of our synthetic target model.

These results confirm that VAEs can model pairwise epistasis in protein
sequence datasets, since they generate pairwise mutational covariances that
are correlated with the target values, even in the absence of explicit
constraints for reproducing these statistics. However, they scale down the
strength of pairwise covariances in both the synthetic and natural analyses
and the correlation with the target is lower than 1. Mi3, in contrast, is
constrained by design to fit the pairwise covariance scores and does so nearly
perfectly.

### Higher order marginal statistics

A more stringent test of GPSM generative capacity is to measure the model’s
ability to reproduce sequence covariation involving more than two positions,
or higher-order covariation. We characterize these higher-order covariation
patterns in the target MSA and GPSM-generated evaluation MSAs by computing the
frequency of non-contiguous amino acid $n$-tuples, or higher-order marginals
(HOMs) corresponding to subsequences, and compare their frequency in each MSA
to corresponding values in the target MSAs. For increasing values of $n$ the
number of possible $n$-tuple combinations increases rapidly, requiring
increasingly large evaluation MSAs to accurately estimate the frequency of
individual $n$-tuples. For this reason, we limit $n$-tuple length to $n\leq
10$ and only compute a limited subset of all possible position sets for each
$n$. For each $n$ we randomly choose 3K position sets, compute the frequencies
of the top twenty most frequent $n$-tuples for each corresponding position set
in the target and evaluation MSAs, as these are well sampled, and for each
position set compute the Pearson correlation between these top twenty
frequencies. We then average the correlation values for each $n$ over all
position sets. We call this metric the Pearson correlation $r_{20}$ [16]. In
this test, estimation error is non-zero because of the extremely large MSAs
required to compute $n$-tuple frequencies, particularly for high $n>5$ (see
Supplementary Information). We can predict the estimation error caused by
finite sampling in the evaluation MSAs by computing the $r_{20}$ scores
between two non-overlapping MSAs generated by the synthetic target model,
which are of the same size as our evaluation MSAs.

In Fig. 2, Bottom Row, we plot the HOM $r_{20}$ for varying $n$. The expected
estimation error (black line) represents a generative capacity upper bound,
giving the highest measurable $r_{20}$ given the evaluation MSA size of 6M for
the synthetic analysis and 10K for the natural analysis. The $r_{20}$ for Mi3
fit to 1M training sequences is very close to the validation upper-bound for
all $n$, suggesting it has accurately fit the synthetic target distribution
and its specification error is close to zero (Fig. 2(d)). This is expected
since the synthetic target model in this test is a Potts model. With 10K
training sequences, Mi3 $r_{20}$ scores are lower than the 1M result for all
$n$ (Fig. 2(e)), which illustrates that Mi3 is affected by out-of-sample error
for typical dataset sizes, as previously described [46]. Indep has much lower
$r_{20}$ scores than Mi3, as expected, since it does not model pairwise
epistasis by design. Indep’s $r_{20}$ scores are similar across all $r_{20}$
experiments, suggesting that it is not strongly affected by out-of-sample nor
estimation error for this metric. This is expected, because its parameters are
optimized for reproducing single-site frequency statistics only, which can be
accurately estimated even from small training MSAs [46]. The $r_{20}$ scores
for DeepSequence are slightly higher than sVAE’s, but both VAEs remain close
to each other for all training datasets and $n$, falling approximately halfway
between Mi3 and Indep, and generally well below Mi3 at higher orders. Here in
the 1M training MSA experiment, the VAEs’ $r_{20}$ decreases to $\sim$0.4 and
$\sim$0.5 at $n=10$ for sVAE and DeepSequence respectively, reflecting higher
specification error at higher orders (Fig. 2(d)). With 10K synthetic training
sequences, the VAEs’ $r_{20}$ decreases further for all $n$ due to the
addition of out-of-sample error (Fig. 2(e)).

Whereas the pairwise covariation correlations represent a preliminary
indication that VAEs capture epistasis but mispredict its strength, the
$r_{20}$ results reinforce this finding and extend it into higher orders.
Unlike Mi3, the VAEs show specification error even when fit to large datasets
from a model which only contains pairwise epistatic interaction terms (Fig.
2(d)). Because higher-order covariation statistics are constrained by the
pairwise statistics, and the VAEs mispredict the pairwise statistics, we
expect that the VAEs will exhibit specification error for higher-order
epistasis, even though our synthetic test does not address this issue
directly. When considered together with the $r_{20}$, the pairwise covariance
correlations reveal a novel insight, which is that when trying to gauge GPSM
generative capacity at higher orders of covariation, the standard pairwise
statistics alone can be misleading. The relative magnitudes of $r_{20}$
between models at $n=2$ are different at higher $n$, and the performance
decrease as $n$ increases is more severe for the VAEs than Mi3 (Fig. 2, Bottom
Row). In the natural analysis, Mi3 performs as close to the target as is
measurable within error given the limited amount of data (Fig. 2(f)). Here,
the small 10K target MSA causes large estimation error, making it more
difficult to distinguish the performance of the different models. We emphasize
that the $r_{20}$ performance decrease for Mi3 when tested against the
synthetic and natural MSA targets can be accounted for entirely by out-of-
sample error and estimation error, consistent with the hypothesis that
pairwise interaction terms are sufficient to model protein sequence variation.

### Hamming distance distributions

We next evaluate the pairwise Hamming distance distribution metric
$d(S,S^{\prime})$. The Hamming distance between two protein sequences is the
number of amino acids that are different between them, and we obtain a
distribution for an MSA by comparing all pairs of sequences. Because it
characterizes the range of sequence diversity in an MSA, recapitulation of the
Hamming distance distribution has been used in the past as a measure for GPSM
performance [2, 16, 36, 37]. In Fig. 3, we compare the pairwise Hamming
distance distribution for each GPSM to that of the target distribution,
computed with evaluation and target MSAs of 10K sequences each. To quantify
the difference between the GPSM and target distributions for this metric, we
use the total variation distance ($\operatorname{TVD}$) [47], which equals 1
when the distributions have no overlap and is 0 when they are identical,
defined by $\operatorname{TVD}[f,g]=1/2\int|f(x)-g(x)|dx$.

All models reproduce the mode Hamming distance of $\sim$179\. For Mi3, we
report the same $\operatorname{TVD}=$0.007 when trained on either 1M (Fig.
3(a)) or 10K (Fig. 3(b)) synthetic sequences, showing negligible specification
error, as expected. When trained on 10K natural sequences, Mi3
$\operatorname{TVD}$ increases to 0.012 (Fig. 3(c)), for reasons discussed
further below. Indep severely underestimates the probability of both low and
high Hamming distances, as observed at the distribution tails, with
$\operatorname{TVD}\sim$0.24 across all experiments. The VAEs perform in
between Mi3 and Indep, but much closer to Mi3 than Indep with respect to
$\operatorname{TVD}$. Performance differences across all GPSMs for this metric
indicate that out-of-sample error has a consistent and detectable, though very
small, effect on the fundamental sequence diversity of artificial GPSM-
generated MSAs. That Mi3 and the VAEs are highly performant and comparable to
each other, but not Indep, corroborates our earlier findings that epistasis is
relevant to accurate modelling of protein sequence diversity (Fig. 2).
However, because Indep performs much closer to Mi3 and the VAEs for this
metric than on any other, and also because this metric cannot discriminate
well between Mi3 and the VAEs, we suspect that reproducing the Hamming
distance distribution is a much easier hurdle for GPSMs than is reproducing
higher-order covariation. This shortcoming of the standard Hamming distance
distribution metric becomes apparent when these results are compared to those
of our novel metric, $r_{20}$, which does show a significant gap in generative
capacity between Mi3 and the VAEs at higher orders (Fig. 2, Bottom Row).

To emphasize the decay of the tails, we rescale all the distributions by their
maxima and re-center them around their modes to give them the same peak, and
then plot them on a log-log scale (Fig. 3, Bottom Row). The relevance of the
distributions’ tails lies in their power-law behavior as they approach zero,
where the function’s exponent is related to the intrinsic dimension of the
dataset and therefore to the number of informative latent factors needed to
explain the data [36, 37, 48]. A well-specified GPSM ought to reproduce this
exponent, and therefore the tail’s decay, since it is a topological property
intrinsic to the dataset and independent from the particular choice of
variables used to describe the probability density [48]. There is a trend of
slightly decreasing generative capacity as training samples decrease, which is
detectable only here in the log-scaled Hamming distribution (Fig. 3, Bottom
Row). In this modified rendering of the Hamming distance distribution,
differences in GPSM generative capacity can be observed at both low (Left
Tail) and high (Right Tail) sequence diversity. The Mi3 distribution closely
overlaps the target distribution with both 1M (Fig. 3(d)) and 10K (Fig. 3(e))
synthetic training sequences. In the 10K natural experiment, Mi3 deviates
noticeably from the target on the left tail (Fig. 3(f)), which represents less
evolutionarily diverged sequences. This could be an artifact of the
phylogenetic relationships between sequences present in the natural dataset,
which may have been incompletely removed by our phylogenetic filtering step
for this dataset (see Methods, Supplementary Information), or it could be due
to estimation error in measuring the target distribution, as only 10K target
sequences are available to estimate the black line in the natural analysis. As
before, Indep performance is consistently low across all experiments. The
VAEs’ performance at low sequence diversity (Fig. 3, Bottom Row, Left Tails)
decreases for smaller training dataset size.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Pairwise Hamming distance distributions. These plots illustrate
whether sequences generated from the GPSMs reproduce the overall sequence
diversity of their respective targets. Mi3 (cyan), DeepSequence (green), sVAE
(red), and Indep (gray) distributions are compared to target distribution
(dotted black). GPSMs were trained on 1M synthetic (a, d), 10K synthetic (b,
e), or 10K natural (c, f) sequences from the corresponding target
distribution. All Hamming distributions were computed from 30K-sequence MSAs,
except for the natural target, which was computed from a 10K-sequence target
MSA due to data limitations. Top Row, Hamming distances $d$ (x-axis) are shown
about the mode, and frequency $f$ is normalized as a fraction of total
(y-axis). Mi3 perfectly matches the target distribution from a to c, whereas
DeepSequence and sVAE overlap each other and share a mode with a slightly
higher frequency than the target and Mi3. Indep has a slightly lower mode than
all the other GPSMs, and a much higher mode frequency. Because generative
capacity for all GPSMs is unchanged from a to c, none of them are sensitive to
training or evaluation sample size for this metric. Bottom Row, Re-scaled
logarithmic Hamming distance distributions better discriminate between GPSMs
with respect to generative capacity than the normal Hamming distance
distribution. Before being log-scaled, the Hamming distances $d$ are
normalized by the mode $d_{Mo}$ (x-axis), and frequencies $f$ are normalized
by the maximum Hamming distance $f_{max}$ (y-axis). This transformation
highlights minute differences between distributions at low frequencies in the
tails of the distributions on the left- and right-hand sides. From d to e,
DeepSequence and sVAE appear sensitive to training sample size for this metric
at the log-log scale.

### Statistical Energy Correlations

A fourth metric we use to evaluate generative capacity is the statistical
energy $E(S)$ of individual sequences in the dataset, which we express using
the negative logarithm of the predicted sequence probability $p(S)$, where
$E(S)=-\log p(S)$. $E(S)$ can be computed analytically for Mi3 and Indep, and
estimated for VAE models by importance sampling (see Methods).

This statistic directly evaluates accuracy of the GPSM distribution values
from $p_{\theta}(S)$ for a limited number of individual sequences, which has
been used to validate GPSMs by comparison to corresponding experimental
fitness values [38, 20, 49, 16]. In Fig. 4, we compare artificial statistical
energies from the GPSM distribution $p_{\theta}(S)$ to those of the target
distribution $p^{0}(S)$ for a 1K test MSA generated from $p^{0}(S)$. This
measurement cannot be performed for the natural 10K experiment because
$p^{0}(S)$ is unknown for the natural data, so we present results only for the
1M and 10K synthetic experiments. We use the 1M (Fig. 4, Left Column) and 10K
(Fig. 4, Right Column) synthetic training MSA sizes, and quantify GPSM
generative capacity for this metric by the Pearson correlation coefficient
$\rho(\\{E(S)\\},\\{\hat{E}(S)\\})$ between synthetic target energies $E(S)$
and GPSM energies $\hat{E}(S)$. Mi3 reproduces the synthetic target
distribution at both training MSA sizes. Because Mi3 should have very low
specification error on the synthetic target, as it is well specified by
design, the small amount of error must be due to remaining out-of-sample or
numerical errors. As expected, Indep poorly reproduces the target values, with
$\rho=0.6$ for both MSA training sizes. The VAEs exhibit slightly larger
specification error than Mi3 on the 1M training set with correlation of
$\rho=0.94$, and exhibit further out-of-sample error on the 10K training set
with $\rho=0.89$.

Juxtaposing the $E(S)$ results to the $r_{20}$ results reiterates the striking
insight of our work, which is that despite the utility of standard metrics for
measuring sequence statistics, they say little about a GPSM’s ability to
capture higher-order covariation, and therefore necessarily higher-order
epistasis. If considered by itself, the $E(S)$ metric would indicate that both
VAEs have generative capacity close to that of Mi3, and even Indep could be
said to have a large amount of generative capacity, in spite of capturing no
covariation by design. But according to $r_{20}$, which directly measures the
higher-order marginals, the VAEs and Indep are comparable to Mi3 only at the
pairwise level. Critically, at higher orders, the performance difference
between the VAEs to each other remains relatively unchanged, but is
significantly lower than Mi3, with Indep lower still. This suggests that
$E(S)$, just as with the Hamming distance metric and pairwise covariance
correlation, represents an easier, and perhaps different, hurdle for GPSMs
than measuring generative capacity by the ability to capture higher-order
sequence covariation, as done uniquely by $r_{20}$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4: Statistical energy correlations. Statistical energies $E(S)$ of 1K
synthetic test sequences from the target distribution as evaluated by Mi3
(cyan), DeepSequence (green), sVAE (red), Indep (gray). Each GPSM was trained
on 1M (a, c, e, g) or 10K (b, d, f, h) sequences from the synthetic target
distribution. For each scatterplot, Pearson correlation coefficient $\rho$ was
computed between each GPSM’s statistical energy and that of the synthetic
target distribution for each sequence. These statistical energies are unit-
less, and so we have shifted the data along the y-axis to compare $E(S)$ on
the same scale. Only Indep is insensitive to decreased training sample size
for this metric.

## Conclusions

In this study we reveal the steep challenges and limits entailed by current
measurements of generative capacity in GPSMs trained on either synthetic or
currently available natural sequence datasets. Recent state-of-the-art GPSM
studies have benchmarked their models by comparing $p_{\theta}(S)$ to
experimental fitness values from deep mutational scans [20, 22], or by
generating artificial sequences that appear to fold into realistic structures
based on in silico folding energy [25]. However, these strategies for model
evaluation present their own challenges. Point mutation fitness experiments in
practice may measure particular contributors to fitness, including replicative
capacity, drug resistance, protein stability, and enzymatic activity [20], and
are subject to significant experimental error and other limitations, e.g.
those imposed by conservation [50]. Likewise, in silico computational
chemistry approaches rely on energy functionals that may lack the precision
needed to meaningfully discriminate between highly similar sequences [51].
Additionally, neither protein function nor fitness rely exclusively on the
thermodynamic stability of static native structures, but also on the protein’s
conformational dynamics [52, 53, 54, 55, 56], which are not fully described by
folding energy values alone [57]. This could mean that despite generating
sequences with realistic in silico folding energy, a GPSM may still not be
capturing crucial higher-order epistatic effects. Neither point mutation
fitness effects, nor in silico folding energy estimations, are directly
related to mutational covariation statistics observed in an MSA in the sense
that they do not check if subsets of covarying amino acids in specific
positions present in the target MSA are indeed present in the GPSM-generated
evaluation MSA. Our novel $r_{20}$ metric uniquely delivers that
functionality, emphasizing higher-order covariation where previous studies
rarely go beyond the pairwise level [2, 16, 58, 34, 35].

Benchmarking coevolution-based protein sequence models in data rich and data
poor regimes, as done here, is an effective method for ascertaining where
data-driven effects stop, and algorithmic failure begins [40]. Due to the
limited availability of natural protein sequence data, this line is inherently
blurred in the natural analysis, as we demonstrate across all four of our
generative capacity metrics. But in our synthetic analysis track, we have
demonstrated the extent to which VAEs, with different implementations, can
capture higher-order covariation at orders between three and ten when the
target distribution is known, its statistical properties are measurable with a
high degree of certainty, and major forms of error are removed, minimized, or
accounted for. When given a large number of training and target sequences, we
found both VAEs’ generative capacity to be between that of a site-independent
model (Indep) and a pairwise Hamiltonian (Mi3) for all measurements. In the
synthetic $r_{20}$ tests, our results show that both VAEs’ generative capacity
is well below Mi3, raising questions about whether VAEs can capture higher-
order epistasis significantly better than a pairwise Potts model. In our
synthetic analysis the target distribution is a Potts model, and therefore we
expect Mi3 to fit the target distribution well by design. However, we find Mi3
also outperforms the VAEs where we do not have this expectation, such as on
the natural target distribution and on a target distribution specified by
sVAE, as shown in Supplementary Information, suggesting Mi3 generally
outperforms VAEs on protein sequence data with respect to generative capacity.

The Hamming distance distributions, pairwise covariation correlations, and
statistical energy correlations are standard metrics that have been used in
the past to measure GPSM accuracy, but we find that they can be inadequate or
misleading indicators of a GPSM’s ability to capture covariation at higher
orders. Taken together, our results suggest that, of the metrics we tested,
only $r_{20}$ provides the granularity needed to discriminate between
different GPSM’s ability to model higher-order epistasis, as it directly tests
the model’s ability to capture higher-order covariation.

Although our results suggest VAEs are less effective for capturing higher-
order epistasis than pairwise Potts models, VAEs have demonstrated utility in
unsupervised learning and clustering. One VAE-GPSM, “BioSeqVAE”, has generated
artificial sequences that share a “hallucinated” homology to natural proteins
in the training set, which could mean that their folded structures would
perform similar functions to their hallucinated natural homologs [23].
Another, “PEVAE”, has shown that a VAE-GPSM’s latent space captures
phylogenetic relationships [21] better than PCA [59] and t-SNE [60]. These
VAE-GPSMs furnished a latent space that immediately allowed for function-based
protein classification, a benefit unavailable to pairwise Potts models without
some effort.

The causes of VAE performance limitations are actively being investigated in
literature from multiple directions. One line of inquiry involves a VAE
phenomenon known as “posterior collapse” (see Supporting information), in
which some dimensions of the VAE latent space become insensitive to the input
data. Studies of this phenomenon have led to some insights into VAE behavior,
for instance that in some situations the VAE likelihood can contain spurious
local maxima[61], and many different heuristic strategies to understand and
avoid this phenomenon have been suggested [62, 63, 64]. In Supporting
Information we test that sVAE, used in our main results, does not exhibit
posterior collapse, though we can trigger it for sVAE architectures with more
than 7 latent dimensions in the bottleneck layer. Another line of inquiry
relates to assumptions typically made about the metric and topology of the
latent space, suggesting that the commonly used Euclidean metric space or
Gaussian prior distribution may not best describe particular datasets, for
instance because of an effect called “manifold mismatch” [65, 66]. Techniques
closely related to the VAE, such as the “WAE”, are also being investigated as
alternatives to the VAE, with different performance characteristics[67]. Among
these competing frameworks for exploring VAE behavior, we have tested VAE
implementations which are currently used to generatively model protein
sequences. With more nuanced latent variable models, and with better
understanding of protein sequence embeddings, perhaps GPSM generative capacity
could extend beyond what has been demonstrated here with state-of-the-art
VAEs.

Our epistasis-oriented methodology focuses on measuring higher-order
covariation, with the potential for broad applicability to various
sequentially ordered data. $r_{20}$-like measurements become possible when the
data are sufficient in number, and the correlation structures between
elements, both within and across samples, are statistically detectable and
meaningful in some context, be it visual, biophysical, or linguistic. The
convergence between data categories such as images, proteins, and language
with respect to generative modelling evaluation offers the exciting
opportunity of a wider, interdisciplinary audience for the work proposed here.
Conversely, further development of sophisticated, data-intensive, and direct
generative capacity metrics of GPSMs could reveal nuances of the correlation
structure of protein sequence datasets that distinguish them from other
datasets, helping to explain why $r_{20}$-like metrics can detect higher-order
covariation, whereas other metrics cannot. Our work represents not only a
revision of currently prevailing paradigms of GPSM benchmarking, but also a
challenge to generative protein sequence modelling more broadly, to consider
how epistasis and direct higher-order covariation metrics like $r_{20}$ can
inform their models and results.

## Methods

### Sequence dataset preparation

An outline of our sequence processing is shown in Fig.1. In Stage 1 we obtain
sequences. For the natural analysis, we use an MSA of the kinase protein
superfamily which we have previously curated using sequences from the
Uniprot/TREMBL database [41]. This MSA is composed of $\sim$20K sequences of
length 232, obtained by filtering a larger set of $\sim$291K sequences to
remove any sequences with more than 50% sequence identity to another (see
Supplementary Information). For the synthetic analysis, we treat the Potts
Hamiltonian model trained on this natural protein kinase MSA as the target
model or distribution, and generate MSAs from it which serve as training and
target MSAs in our synthetic tests.

In Stage 2 we split the sequences into training and target MSAs. For the
natural analysis, we randomly divide the 20K natural sequences into 10K
training and 10K target MSAs. For the synthetic analysis we generate both
training MSAs used to train the synthetic models, and target MSAs used as
reference, and ensure that the synthetic training and target sequences are
non-overlapping sets, even though they come from the same target model.

In Stage 3, the processed MSAs are then one-hot encoded and used to train the
the GPSMs. In Stage 4, we generate evaluation MSAs from each model, and in
Stage 5 we perform our generative capacity measurements by comparing the
artificial MSAs to the appropriate target MSA.

### Mi3

The Mi3 model is a pairwise Potts Hamiltonian model fit to sequence data using
the “Mi3-GPU” software we have developed previously [33], which performs
“inverse Ising inference” to infer parameters of Potts models using a Markov-
Chain Monte-Carlo (MCMC) algorithm which entails very few approximations. This
software allows us to fit statistically accurate Potts models to MSA data. We
have examined Mi3’s generative capacity and out-of-sample error in earlier
work [33, 46], which we summarize here.

A Potts model is the maximum entropy model for $p(S)$ constrained to reproduce
the bivariate marginals $f^{ij}_{\alpha\beta}$ of an MSA, i.e. the frequency
of amino acid combination $\alpha,\beta$ ,at positions $i,j$. The probability
distribution $p_{\theta}(S)$ for the Potts model takes the form

$p_{\theta}(S)=\frac{e^{-E(S)}}{Z}\quad\text{with}\quad
E(S)=\sum_{i}^{L}h^{i}_{s_{i}}+\sum_{i<j}J^{ij}_{s_{i}s_{j}}$ (1)

where $Z$ is a normalization constant, $Z=\sum_{S}e^{-E(S)}$, and “coupling”
$J^{ij}_{\alpha\beta}$ and “field” $h^{i}_{\alpha}$ parameters are compactly
denoted by the vector $\theta=\\{h^{i}_{\alpha},J^{ij}_{\alpha\beta}\\}$. The
number of free parameters of the model (non-independent couplings and fields)
is equal to the number of non-independent bivariate marginals, which can be
shown to be $\frac{L(L-1)}{2}(q-1)^{2}+L(q-1)$ for $q$ amino acids [14], which
is $\sim$10.7M parameters for our model. This implies that the Potts model is
well specified to reproduce the bivariate marginals when generating sequences
from $p_{\theta}(S)$.

The Mi3 model inference procedure maximizes the log-likelihood with
regularization. Maximizing the Potts log-likelihood can be shown to be
equivalent to minimizing the difference between the dataset MSA bivariate
marginals $f^{ij}_{\alpha\beta}$ and the model bivariate marginals of
sequences generated from $p(S)$. To account for finite sampling error in the
estimate of $f^{ij}_{\alpha\beta}$ for an MSA of $N$ sequences we add a small
pseudocount of size $1/N$, as described previously [46]. We also add a
regularization penalty to the likelihood affecting the coupling parameters
$J^{ij}_{\alpha\beta}$ to bias them towards 0, of form
$\lambda\sum\operatorname{SCAD}(J^{ij}_{\alpha\beta},\lambda,\alpha)$ using
the SCAD function which behaves like $\lambda|J^{ij}_{\alpha\beta}|$ for small
$J^{ij}_{\alpha\beta}$ but gives no bias for large $J^{ij}_{\alpha\beta}$[68],
using a small regularization strength of $\lambda=0.001$ for all inferences
which causes little model bias.

To generate synthetic MSAs from the Potts model we use MCMC over the trial
distribution $p_{\theta}(S)$ until the Markov-Chains reach equilibrium, as
described previously [33]. We can directly evaluate $E(S)$ as the negative
log-probability of any sequences for the Mi3 model using equation 1 up to a
constant $Z$, and this constant can be dropped without affecting our results.

### Indep

The Indep model is the maximum entropy model for $p(S)$ constrained to
reproduce the univariate marginals of an MSA, and is commonly called a “site-
independent” model because the sequence variations at each site are
independent of the variation at other sites. Because it does not fit the
bivariate marginals, it cannot model covariation between positions. It takes
the form

$p_{\theta}(S)=\frac{e^{-E(S)}}{Z}\quad\text{with}\quad
E(S)=\sum_{i}^{L}h^{i}_{s_{i}}$ (2)

where $Z$ is a normalization constant, $Z=\sum_{S}e^{-E(S)}$, and “field”
parameters $h^{i}_{\alpha}$ for all positions $i$ and amino acids $\alpha$ are
compactly referred to by the vector $\theta=\\{h^{i}_{\alpha}\\}$. The fields
of the Indep model generally have different values from the fields of the
Potts model. Unlike for the Potts model, maximum likelihood parameters can be
determined analytically to be $h^{i}_{\alpha}=-\log f^{i}_{\alpha}$ where
$f^{i}_{\alpha}$ are the univariate marginals of the dataset MSA. When fitting
the Indep model to a dataset MSA of $N$ sequences, we add a pseudocount of
$1/N$ to the univariate marginals to give model marginals
$\hat{f}^{i}_{\alpha}$ , to account for finite sample error in the univariate
marginal estimates. The model distribution simplifies to a product over
positions, as $p_{\theta}(S)=\sum\hat{f}^{i}_{s_{i}}$. The number of
independent field parameters is $L(q-1)$ which equals 4.6K parameters for our
model.

To generate sequences from the independent model we independently generate the
residues at each position $i$ by a weighted random sample from the marginals
$f^{i}_{\alpha}$, and we directly evaluate the log probability of each
sequence $E(S)$ from equation 2.

### VAEs

The standard variational autoencoder (sVAE) is a deep, symmetrical, and
undercomplete autoencoder neural network composed of a separate encoder
$q_{\phi}(Z|S)$ and decoder $p_{\theta}(S|Z)$ [69], which map input sequences
$S$ to regions within a low-dimensional latent space $Z$ and back. The
probability distribution for the sVAE is defined as

$p_{\theta}(S)=\int p_{\theta}(S|Z)p(Z)dZ$ (3)

where the latent space distribution is a unit Normal distribution,
$p(Z)=\mathcal{N}[0,1](Z)$. Training of a VAE can be understood as
maximization of the dataset log-likelihood with the addition of a Kullback-
Leibler regularization term
$\operatorname{D_{KL}}[q_{\phi}(Z|S),p_{\theta}(Z|S)]$, where
$p_{\theta}(Z|S)$ is the posterior of the decoder [29, 30].

sVAE’s architecture is “vanilla” [70], meaning it is implemented in a standard
way and its behavior is meant to be representative of VAEs generally. Nearly
identical sVAE architectures have been used as VAE-GPSMs in recent studies
[22]. sVAE’s encoder and decoder are implemented without advancements such as
convolutional layers [20], multi-stage training [71], disentanglement learning
[70], Riemannian Brownian motion priors [72], and more. This allows us to
directly interrogate the assumptions and performance of standard variational
autoencoding with respect to the training and evaluation of GPSMs in this
work.

Our encoder and decoder have 3 layers each and employ standard normalization
and regularization strategies [73, 74]. sVAE’s latent bottleneck layer has 7
nodes, and the model in total has 2.7M inferred parameters. The input layer of
the encoder accepts a one-hot encoded sequence and the decoder’s output layer
values can be interpreted as a Bernoulli distribution of the same dimensions
as a one-hot encoded sequence. We have tested various sVAE architectures and
hyperparameters with our datasets, as well as DeepSequence as described below,
and found qualitatively similar generative capacity results.

To generate a sequence from sVAE, we generate a random sample in latent space
from the latent distribution $p(Z)$ pass this value to the decoder to obtain a
Bernoulli distribution, from which we sample once. To evaluate the negative
log-probability of a sequence $E(S)$ we use importance sampling, averaging
over 1K samples from the latent distribution $q_{\phi}(Z|S)$ [21]. Other
publications use the Evidence Lower Bound (ELBO) estimate as an approximation
of the negative log-probability [20], and we have verified that the ELBO and
the negative log-probability are nearly identical in our tests and have equal
computational complexity.

In addition to sVAE, we test the DeepSequence VAE [20]. We use the default
inference parameters, and use the “SVI” inference implementation which uses a
“variational Bayes” inference technique as an extension of the sVAE inference
method. Because DeepSequence is designed to output the ELBO rather than the
negative log-probability for each sequence, we use the ELBO as an
approximation of $E(S)$, and estimate it using an average over 1K samples. To
generate sequences we use the same strategy as for sVAE.

## References

  * [1] Levy, R. M., Haldane, A. & Flynn, W. F.  Potts hamiltonian models of protein co-variation, free energy landscapes, and evolutionary fitness.  _Current Opinion in Structural Biology_ 43, 55–62 (2017).
  * [2] Cocco, S., Feinauer, C., Figliuzzi, M., Monasson, R. & Weigt, M.  Inverse statistical physics of protein sequences: a key issues review.  _Reports on Progress in Physics_ 81, 032601 (2018).  Publisher: IOP Publishing.
  * [3] Tubiana, J., Cocco, S. & Monasson, R.  Learning compositional representations of interacting systems with restricted boltzmann machines: Comparative study of lattice proteins.  _Neural Computation_ 31, 1671–1717 (2019).  URL https://doi.org/10.1162/neco_a_01210.
  * [4] Lapedes, A. S., Giraud, B., Liu, L. & Stormo, G. D.  _Correlated mutations in models of protein sequences: phylogenetic and structural effects_ (Institute of Mathematical Statistics, 1999).  Pages: 236-256 Publication Title: Statistics in molecular biology and genetics.
  * [5] Weigt, M., White, R. A., Szurmant, H., Hoch, J. A. & Hwa, T.  Identification of direct residue contacts in protein–protein interaction by message passing.  _Proceedings of the National Academy of Sciences_ 106, 67–72 (2009).  Publisher: National Academy of Sciences Section: Physical Sciences.
  * [6] Haldane, A., Flynn, W. F., He, P., Vijayan, R. S. K. & Levy, R. M.  Structural propensities of kinase family proteins from a Potts model of residue co-variation.  _Protein Science: A Publication of the Protein Society_ 25, 1378–1384 (2016).
  * [7] Domingo, J., Baeza-Centurion, P. & Lehner, B.  The Causes and Consequences of Genetic Interactions (Epistasis).  _Annual Review of Genomics and Human Genetics_ 20, 433–460 (2019).  Publisher: Annual Reviews.
  * [8] Noel, J. K., Morcos, F. & Onuchic, J. N.  Sequence co-evolutionary information is a natural partner to minimally-frustrated models of biomolecular dynamics.  _F1000Research_ 5 (2016).  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4755392/.
  * [9] Morcos, F., Jana, B., Hwa, T. & Onuchic, J. N.  Coevolutionary signals across protein lineages help capture multiple protein conformations.  _Proceedings of the National Academy of Sciences_ 110, 20533–20538 (2013).  URL https://www.pnas.org/content/110/51/20533.  Publisher: National Academy of Sciences Section: Biological Sciences.
  * [10] Sułkowska, J. I., Morcos, F., Weigt, M., Hwa, T. & Onuchic, J. N.  Genomics-aided structure prediction.  _Proceedings of the National Academy of Sciences_ 109, 10340–10345 (2012).  URL https://www.pnas.org/content/109/26/10340.  Publisher: National Academy of Sciences Section: Biological Sciences.
  * [11] Hopf, T. A. _et al._ Mutation effects predicted from sequence co-variation.  _Nature biotechnology_ 35, 128–135 (2017).
  * [12] Biswas, A., Haldane, A., Arnold, E. & Levy, R. M.  Epistasis and entrenchment of drug resistance in HIV-1 subtype B.  _eLife_ 8, e50524 (2019).  Publisher: eLife Sciences Publications, Ltd.
  * [13] Schneidman, E., Still, S., Berry, M. J. & Bialek, W.  Network information and connected correlations.  _Phys. Rev. Lett._ 91, 238701 (2003).  URL https://link.aps.org/doi/10.1103/PhysRevLett.91.238701.
  * [14] Bialek, W. & Ranganathan, R.  Rediscovering the power of pairwise interactions.  _arXiv:0712.4397 [q-bio]_ (2007).  ArXiv: 0712.4397.
  * [15] Figliuzzi, M., Barrat-Charlaix, P. & Weigt, M.  How Pairwise Coevolutionary Models Capture the Collective Residue Variability in Proteins?  _Molecular Biology and Evolution_ 35, 1018–1027 (2018).  URL https://doi.org/10.1093/molbev/msy007.  https://academic.oup.com/mbe/article-pdf/35/4/1018/24597926/msy007.pdf.
  * [16] Haldane, A., Flynn, W. F., He, P. & Levy, R. M.  Coevolutionary Landscape of Kinase Family Proteins: Sequence Probabilities and Functional Motifs.  _Biophysical Journal_ 114, 21–31 (2018).
  * [17] Weinreich, D. M., Lan, Y., Wylie, C. S. & Heckendorn, R. B.  Should evolutionary geneticists worry about higher-order epistasis?  _Current opinion in genetics & development_ 23, 700–707 (2013).
  * [18] Haq, O., Andrec, M., Morozov, A. V. & Levy, R. M.  Correlated electrostatic mutations provide a reservoir of stability in hiv protease.  _PLoS Comput Biol_ 8, e1002675EP (2012).
  * [19] Haq, O., Levy, R., Morozov, A. & Andrec, M.  Pairwise and higher-order correlations among drug-resistance mutations in hiv-1 subtype b protease.  _BMC Bioinformatics_ 10, S10 (2009).
  * [20] Riesselman, A. J., Ingraham, J. B. & Marks, D. S.  Deep generative models of genetic variation capture the effects of mutations.  _Nature Methods_ 15, 816–822 (2018).
  * [21] Ding, X., Zou, Z. & Brooks Iii, C. L.  Deciphering protein evolution and fitness landscapes with latent space models.  _Nature Communications_ 10, 5644 (2019).  Number: 1 Publisher: Nature Publishing Group.
  * [22] Sinai, S., Kelsic, E., Church, G. M. & Nowak, M. A.  Variational auto-encoding of protein sequences.  _NeurIPS 2017 MLCB workshop_ (2017).  ArXiv: 1712.03346.
  * [23] Costello, Z. & Martin, H. G.  How to Hallucinate Functional Proteins.  _arXiv:1903.00458 [q-bio]_ (2019).  ArXiv: 1903.00458.
  * [24] Gupta, A. & Zou, J.  Feedback GAN for DNA optimizes protein functions.  _Nature Machine Intelligence_ 1, 105–111 (2019).  Number: 2 Publisher: Nature Publishing Group.
  * [25] Madani, A. _et al._ ProGen: Language Modeling for Protein Generation.  preprint, Synthetic Biology (2020).  URL http://biorxiv.org/lookup/doi/10.1101/2020.03.07.982272.
  * [26] Vig, J. _et al._ {BERT}ology meets biology: Interpreting attention in protein language models.  In _International Conference on Learning Representations_ (2021).  URL https://openreview.net/forum?id=YWtLZvLmud7.
  * [27] Elnaggar, A. _et al._ ProtTrans: Towards Cracking the Language of Life’s Code Through Self-Supervised Deep Learning and High Performance Computing.  _bioRxiv_ 2020.07.12.199554 (2020).  Publisher: Cold Spring Harbor Laboratory Section: New Results.
  * [28] Choromanski, K. M. _et al._ Rethinking Attention with Performers.  In _accepted to ICLR 2021 (oral presentation)_ (2021).  URL https://arxiv.org/abs/2009.14794.
  * [29] Kingma, D. P. & Welling, M.  Auto-Encoding Variational Bayes.  _arXiv:1312.6114 [cs, stat]_ (2014).  URL http://arxiv.org/abs/1312.6114.  ArXiv: 1312.6114.
  * [30] Rezende, D. J., Mohamed, S. & Wierstra, D.  Stochastic Backpropagation and Approximate Inference in Deep Generative Models.  In _International Conference on Machine Learning_ , 1278–1286 (PMLR, 2014).  URL http://proceedings.mlr.press/v32/rezende14.html.  ISSN: 1938-7228.
  * [31] Luce, R. D., Bush, R. R. & Galanter, E. (eds.) _Handbook of mathematical psychology: I._ Handbook of mathematical psychology: I. (John Wiley, Oxford, England, 1963).  Pages: xiii, 491.
  * [32] Frawley, W. J.  _International Encyclopedia of Linguistics_ (Oxford University Press, 2003).  Publication Title: International Encyclopedia of Linguistics.
  * [33] Haldane, A. & Levy, R. M.  Mi3-GPU: MCMC-based inverse Ising inference on GPUs for protein covariation analysis.  _Computer Physics Communications_ 107312 (2020).
  * [34] Shimagaki, K. & Weigt, M.  Selection of sequence motifs and generative Hopfield-Potts models for protein families.  _Physical Review E_ 100, 032128 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevE.100.032128.  Publisher: American Physical Society.
  * [35] Hawkins-Hooker, A. _et al._ Generating functional protein variants with variational autoencoders.  _PLOS Computational Biology_ 17, e1008736 (2021).  URL https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008736.  Publisher: Public Library of Science.
  * [36] Facco, E., Pagnani, A., Russo, E. T. & Laio, A.  The intrinsic dimension of protein sequence evolution.  _PLOS Computational Biology_ 15, e1006767 (2019).  Publisher: Public Library of Science.
  * [37] Granata, D. & Carnevale, V.  Accurate Estimation of the Intrinsic Dimension Using Graph Distances: Unraveling the Geometric Complexity of Datasets.  _Scientific Reports_ 6, 31377 (2016).
  * [38] Riesselman, A. _et al._ Accelerating Protein Design Using Autoregressive Generative Models.  _bioRxiv_ 757252 (2019).  Publisher: Cold Spring Harbor Laboratory Section: New Results.
  * [39] El-Gebali, S. _et al._ The Pfam protein families database in 2019.  _Nucleic Acids Research_ 47, D427–D432 (2019).  Publisher: Oxford Academic.
  * [40] AlQuraishi, M.  ProteinNet: a standardized data set for machine learning of protein structure.  _BMC Bioinformatics_ 20, 311 (2019).  URL https://doi.org/10.1186/s12859-019-2932-0.
  * [41] Consortium, T. U.  UniProt: a worldwide hub of protein knowledge.  _Nucleic Acids Research_ 47, D506–D515 (2019).  URL https://academic.oup.com/nar/article/47/D1/D506/5160987.  Publisher: Oxford Academic.
  * [42] Burnham, K. P. & Anderson, D. R.  _Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach_ (Springer-Verlag, New York, 2002), 2 edn.
  * [43] Yaser S. Abu-Mostafa, Malik Magdon-Ismail, Hsuan-Tien Lin-Learning From Data_ A short course-AMLBook.com (2012).pdf | Statistical Classification | Machine Learning.
  * [44] Everitt, B. S. & Skrondal, A.  The Cambridge Dictionary of Statistics.  _Cambridge University Press_ 480 (2010).
  * [45] Mohri, M., Rostamizadeh, A. & Talwalkar, A.  _Foundations of Machine Learning, second edition_ (MIT Press, 2018).  Google-Books-ID: dWB9DwAAQBAJ.
  * [46] Haldane, A. & Levy, R. M.  Influence of multiple-sequence-alignment depth on Potts statistical models of protein covariation.  _Physical Review. E_ 99, 032405 (2019).
  * [47] Levin, D. A. & Peres, Y.  _Markov Chains and Mixing Times_ (American Mathematical Soc., 2017).  Google-Books-ID: f208DwAAQBAJ.
  * [48] Ansuini, A., Laio, A., Macke, J. H. & Zoccolan, D.  Intrinsic dimension of data representations in deep neural networks.  In _Advances in Neural Information Processing Systems_ , 6111–6122 (2019).
  * [49] Figliuzzi, M., Jacquier, H., Schug, A., Tenaillon, O. & Weigt, M.  Coevolutionary Landscape Inference and the Context-Dependence of Mutations in Beta-Lactamase TEM-1.  _Molecular Biology and Evolution_ 33, 268–280 (2016).
  * [50] Haddox, H. K., Dingens, A. S., Hilton, S. K., Overbaugh, J. & Bloom, J. D.  Mapping mutational effects along the evolutionary landscape of HIV envelope.  _eLife_ 7, e34420 (2018).  Publisher: eLife Sciences Publications, Ltd.
  * [51] Alley, E. C., Khimulya, G., Biswas, S., AlQuraishi, M. & Church, G. M.  Unified rational protein engineering with sequence-based deep representation learning.  _Nature Methods_ 16, 1315–1322 (2019).  URL https://www.nature.com/articles/s41592-019-0598-1.  Number: 12 Publisher: Nature Publishing Group.
  * [52] Wei, G., Xi, W., Nussinov, R. & Ma, B.  Protein Ensembles: How Does Nature Harness Thermodynamic Fluctuations for Life? The Diverse Functional Roles of Conformational Ensembles in the Cell.  _Chemical Reviews_ 116, 6516–6551 (2016).  URL https://doi.org/10.1021/acs.chemrev.5b00562.  Publisher: American Chemical Society.
  * [53] Anishchenko, I., Ovchinnikov, S., Kamisetty, H. & Baker, D.  Origins of coevolution between residues distant in protein 3D structures.  _Proceedings of the National Academy of Sciences_ 114, 9122–9127 (2017).  URL https://www.pnas.org/content/114/34/9122.  Publisher: National Academy of Sciences Section: Biological Sciences.
  * [54] Sailer, Z. R. & Harms, M. J.  Molecular ensembles make evolution unpredictable.  _Proceedings of the National Academy of Sciences_ 114, 11938–11943 (2017).  URL https://www.pnas.org/content/114/45/11938.  Publisher: National Academy of Sciences Section: Biological Sciences.
  * [55] Petrović, D., Risso, V. A., Kamerlin, S. C. L. & Sanchez-Ruiz, J. M.  Conformational dynamics and enzyme evolution.  _Journal of The Royal Society Interface_ 15, 20180330 (2018).  URL https://royalsocietypublishing.org/doi/full/10.1098/rsif.2018.0330.  Publisher: Royal Society.
  * [56] Nussinov, R., Tsai, C.-J. & Jang, H.  Protein ensembles link genotype to phenotype.  _PLOS Computational Biology_ 15, e1006648 (2019).  URL https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006648.  Publisher: Public Library of Science.
  * [57] Campitelli, P., Modi, T., Kumar, S. & Ozkan, S. B.  The Role of Conformational Dynamics and Allostery in Modulating Protein Evolution.  _Annual Review of Biophysics_ 49, 267–288 (2020).  Publisher: Annual Reviews.
  * [58] Haq, O., Levy, R. M., Morozov, A. V. & Andrec, M.  Pairwise and higher-order correlations among drug-resistance mutations in HIV-1 subtype B protease.  _BMC Bioinformatics_ 10, S10 (2009).
  * [59] Bishop, C. M.  Pattern recognition and machine learning (2006).  ISBN: 9781493938438 9780387310732 Library Catalog: cds.cern.ch Publisher: Springer.
  * [60] Maaten, L. v. d. & Hinton, G.  Visualizing Data using t-SNE.  _Journal of Machine Learning Research_ 9, 2579–2605 (2008).
  * [61] Lucas, J., Tucker, G., Grosse, R. & Norouzi, M.  Understanding posterior collapse in generative latent variable models.  _ICLR 2019 Workshop DeepGenStruct_ (2019).
  * [62] Dai, B., Wang, Z. & Wipf, D.  The usual suspects? Reassessing blame for VAE posterior collapse.  In III, H. D. & Singh, A. (eds.) _Proceedings of the 37th International Conference on Machine Learning_ , vol. 119 of _Proceedings of Machine Learning Research_ , 2313–2322 (PMLR, 2020).  URL http://proceedings.mlr.press/v119/dai20c.html.
  * [63] Fu, H. _et al._ Cyclical Annealing Schedule: A Simple Approach to Mitigating KL Vanishing.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 240–250 (Association for Computational Linguistics, Minneapolis, Minnesota, 2019).  URL https://www.aclweb.org/anthology/N19-1021.
  * [64] Takahashi, H., Iwata, T., Yamanaka, Y., Yamada, M. & Yagi, S.  Variational Autoencoder with Implicit Optimal Priors.  _Proceedings of the AAAI Conference on Artificial Intelligence_ 33, 5066–5073 (2019).  URL https://ojs.aaai.org/index.php/AAAI/article/view/4439.  Number: 01.
  * [65] Davidson, T. R., Falorsi, L., De Cao, N., Kipf, T. & Tomczak, J. M.  Hyperspherical variational auto-encoders.  _34th Conference on Uncertainty in Artificial Intelligence (UAI-18)_ (2018).
  * [66] Falorsi, L. _et al._ Explorations in Homeomorphic Variational Auto-Encoding.  _ICML 2018 workshop on Theoretical Foundations and Applications of Deep Generative Models_ (2018).  URL http://arxiv.org/abs/1807.04689.  ArXiv: 1807.04689.
  * [67] Tolstikhin, I., Bousquet, O., Gelly, S. & Schoelkopf, B.  Wasserstein auto-encoders.  In _International Conference on Learning Representations_ (2018).  URL https://openreview.net/forum?id=HkL7n1-0b.
  * [68] Fan, J. & Li, R.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American Statistical Association_ 96, 1348–1360 (2001).  URL https://doi.org/10.1198/016214501753382273.  https://doi.org/10.1198/016214501753382273.
  * [69] Charte, D., Charte, F., García, S., del Jesus, M. J. & Herrera, F.  A practical tutorial on autoencoders for nonlinear feature fusion: Taxonomy, models, software and guidelines.  _Information Fusion_ 44, 78–96 (2018).
  * [70] Ding, Z. _et al._ Guided Variational Autoencoder for Disentanglement Learning.  In _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 7917–7926 (2020).  ISSN: 2575-7075.
  * [71] Dai, B. & Wipf, D.  Diagnosing and Enhancing VAE Models.  _arXiv:1903.05789 [cs, stat]_ (2019).  URL http://arxiv.org/abs/1903.05789.  ArXiv: 1903.05789.
  * [72] Kalatzis, D., Eklund, D., Arvanitidis, G. & Hauberg, S.  Variational Autoencoders with Riemannian Brownian Motion Priors.  In _ICML_ (2020).
  * [73] Ioffe, S. & Szegedy, C.  Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift.  In _International Conference on Machine Learning_ , 448–456 (PMLR, 2015).  URL http://proceedings.mlr.press/v37/ioffe15.html.  ISSN: 1938-7228.
  * [74] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I. & Salakhutdinov, R.  Dropout: a simple way to prevent neural networks from overfitting.  _The Journal of Machine Learning Research_ 15, 1929–1958 (2014).

## Acknowledgements

This research was supported by the National Science Foundation, grant number
NSF-GCR 1934848, and the National Institute of Health, grant number
R35-GM132090. This research includes calculations carried out on HPC resources
supported in part by the National Science Foundation through major research
instrumentation grant number 1625061 and by the US Army Research Laboratory
under contract number W911NF-16-2-0189.

## Author contributions statement

F.M., R.M.L, V.C., A.H. conceived the experiments. F.M., V.C., A.H. performed
the experiments. F.M., R.M.L, V.C., A.H. analyzed the results. F.M., V.C.,
A.H. wrote the bulk of the codebase, Q.N. made a contribution to the codebase.
F.M., R.M.L, V.C., A.H. wrote the paper. All authors reviewed the manuscript.