# Few-shot protein generation

Soumya Ram
Massachusetts Institute of Technology
Department of Biological Engineering
Cambridge,MA,USA
Tristan Bepler
Simons Machine Learning Center
New York Structural Biology Center
New York, NY, USA

###### Abstract

We present the MSA-to-protein transformer, a generative model of protein
sequences conditioned on protein families represented by multiple sequence
alignments (MSAs). Unlike existing approaches to learning generative models of
protein families, the MSA-to-protein transformer conditions sequence
generation directly on a learned encoding of the multiple sequence alignment,
circumventing the need for fitting dedicated family models. By training on a
large set of well curated multiple sequence alignments in Pfam, our MSA-to-
protein transformer generalizes well to protein families not observed during
training and outperforms conventional family modeling approaches, especially
when MSAs are small. Our generative approach accurately models epistasis and
indels and allows for exact inference and efficient sampling unlike other
approaches. We demonstrate the protein sequence modeling capabilities of our
MSA-to-protein transformer and compare it with alternative sequence modeling
approaches in comprehensive benchmark experiments.

##  1 Introduction

Protein engineering is the task of mutating proteins in order to achieve a
desired function, and has numerous applications in medicine and sustainability
[1, 2]. However, designing these mutations is challenging, because it is
difficult to infer their functional impact and the search space of possible
sequence variants is combinatorialy large. A given mutation could cause a
disproportionate effect due to being positioned within an active site or long-
range interactions with other amino acids [3]. In addition, introducing
multiple mutations simultaneously can have complex non-linear effects, called
epistasis [4].

Supervised data on the functional impact of protein mutants is also limited.
Acquiring supervised data means performing complicated and costly deep
mutational scanning experiments. The number of variants that can be measured
is limited by the assay throughput. For typical functional activity assays,
this can restrict the number of feasibly measurable variants to hundreds or
less. Generally, these experiments characterize the functional impact of point
mutations, so experimental data on the impact of higher order mutants is even
more rare [5, 6].

In contrast, sequence data is plentiful for natural proteins. The number of
known natural protein sequences has nearly tripled in the last 5 years, and
continues to grow rapidly due to the falling cost of DNA sequencing. This
proliferation has enabled new evolutionary sequenced modeling approaches for
variant function prediction that circumvent the need for functional assay
data. Classically, multiple sequence alignments (MSAs) of related protein
sequences have proven useful for predicting protein contacts [7, 8, 9, 10],
inferring functional effects of mutants [11], and searching databases for
related sequences [12, 13]. Simple generative models learned from MSAs, such
as position-specific scoring matrices (PSSMs) or profile Hidden Markov Models
(pHMMs), are widely used protein sequence models [13].

Recently, deep language models have achieved great success modeling proteins.
Protein embedding models produce embeddings that surpass hand-crafted features
in generalizing to downstream tasks [14, 15, 16, 17, 18]. Family-specific deep
generative sequence models have also demonstrated impressive unsupervised
variant function prediction on deep mutational scanning datasets [6, 19] and
for predicting antibody escape mutations [20]. However, these models must be
trained on a family-by-family basis. [21] recently observed that a transformer
pre-trained as a masked language model on general MSAs is able to predict the
effects of variants without retraining by feeding the family MSA directly as
input to the model. However, this approach assumes an additive model of
multiple variants and is unable to model insertions, limiting its usability.
Inefficient inference in non-autoregressive models also creates challenges for
their practical application to sequence generation.

Here, we introduce a new approach to learning generative models of proteins
based on sequence-to-sequence learning. We newly formulate sequence modeling
as a few-shot learning problem and train a single encoder-decoder model to
accept a protein family represented as an MSAs as input and decode that MSA
into a distribution over sequences from that family. We train the model on
tens of thousands of multiple sequence alignments representing known protein
families and evaluate it on unseen families held out from training. We achieve
state-of-art performance on deep mutational scanning datasets and are able to
model insertions and deletions in addition to substitutions. Our model also
makes no additivity assumptions, allowing it to generalize well to highly
diverse mutants. Unlike previous work, our model offers exact sampling and
inference and crucially, does not require retraining for new families.

##  2 Related Work

Generative sequence models are widely used to represent protein families.
PSSMs model each column in the MSA as independent distributions over amino
acids. Profile HMMs model each amino acid as being generated conditioned on a
hidden state corresponding to each column in the MSA, but this alignment is
considered unobserved when calculating the probability of a new sequence.
These simple models are widely used, because they need to be inferred from
relatively small sets of sequences (often only 10s or 100s) and parameter
inference needs to be performed for each set of proteins of interest. However,
they are unable to model dependencies between amino acids.

Deep generative models have recently been explored for conditional and
unconditional sequence generation, because they can model diverse protein
sequences and capture high-order dependencies between amino acids. [22]
present an autoregressive transformer language model conditioned on taxonomic
and functional labels. [23] develop a graph-based sequence generative model
conditioned on protein tertiary structure. Other work has considered
generative models for adaptively sampling from posterior distributions defined
by function prediction methods, such as [24].

Two of the proposed generative models for proteins require training a VAE on
thousands of sequences for each new family [6, 25]. The third generative model
in the literature evaluated on variant function prediction, MSA Transformer,
also can do few-shot prediction [21]. Compared to MSA Transformer, our
formulation utilizes a decoder that attends over the multiple sequence
alignment to explicitly model sequence probabilities autoregressively. In
contrast, the MSA Transformer is a masked language model which cannot give a
closed-form sequence probability, is unable to generalize to sequences with a
different number of positions from the MSA, and does not permit efficient
sampling.

The addition of a decoder allows for many benefits in inference and sampling.
Because of the decoder, our model is able to produce exact likelihoods for
higher order mutants. Masked zero-shot language models such as MSATransformer
and ESM-1v evaluate likelihoods for higher order mutants by assuming
additivity [21]. However, these proteins tend to have strong epistatic effects
so such an assumption is not ideal [4].

MSA2Prot is also able to generate sequences with insertions as well, allowing
it to model the entire protein landscape. This is especially valuable as
indels comprise of 15-21 percent of human polymorphisms, and inferring their
functional impact is important for understanding genes related to disease
[25].

Sampling-wise, our model produces exact samples while zero-shot masked
language models [21] require Gibbs sampling. Gibbs sampling struggles with
transitioning between high probability states that are connected by low
probability paths. This is the case with many higher order mutants. For
example, double mutants can be beneficial when each of the single mutants are
not because of epistatic interactions. In addition, this issue is compounded
by the fact that directed evolution experiments can require a large number
(possibly thousands) of initial protein samples. Producing this many
independent samples through Gibbs creates a large computational burden.

Thus, MSA2Prot’s formulation models epistatic effects from higher order
mutants and indels, while also allowing for efficient and exact sampling and
evaluation.

Lastly, MSA2Prot can be used for adaptive sampling in settings with scarce
protein fitness measurements as it generalizes from the input MSA without
needing retraining. This is in contrast to methods such as [26] which
continually require the generator to be retrained. Thus, MSA2Prot is able to
combine the data efficiency of simple models such as HMMs with the expressive
generative capability of deep learning models.

##  3 Methods

Let $x_{i,j}$ , $i\in 1...N$ , $j\in 1...M$, denote the amino acid or gap
token at the $j$’th position of the $i$’th sequence in a multiple sequence
alignment (the query MSA) with N sequences and M columns. Let $y_{i}$ , $i\in
1...L$ , be the token for the $i$’th amino acid in the target protein
sequence. Our goal is to learn the probability distribution over target
protein sequences ($y$’s) from protein families represented by their members
in a multiple sequence alignment ($X$’s), $p(y|X)$, given a training set of
MSA, target protein pairs ($X^{k}$,$y^{k}$) where $k\in 1...K$ denotes the
$k$’th pair. Here, target protein $y^{k}$ is a member of the same family as
the sequences in $X^{k}$ that is not itself in $X^{k}$. In order to generalize
across input MSAs, we propose an MSA-encoder, sequence-decoder neural network
architecture. The protein MSA is encoded using an MSA transformer with axial
self-attention on the rows and columns, into context-aware vector
representations, $z_{i,j}$, for each token in the MSA, and then the target
sequence is decoded by attending to the learned MSA representations in a
decoder transformer with cross attention to the encoded MSA.This is
illustrated in Figure 1

Figure 1: MSA2Prot Architecture

###  3.1 The Transformer MSA Encoder

The MSA encoder accepts a multiple sequence alignment as input and returns a
vector representation for each position in the multiple sequence alignment,
$z_{i,j}\in R^{d}$, where $d$ is the dimension of the learned embedding. This
encoder is parameterized as a stack of transformer layers with axial attention
over the rows and columns of the MSA. Before being processed by the
transformer stack, the input tokens are embedded into vectors in $R^{d}$ and
augmented with a random Fourier projection of the column index as a positional
embedding, also in $R^{d}$. No positional embedding is used for the rows of
the MSA (the sequence index), because the ordering of sequences in an MSA is
arbitrary and we wish for the MSA encoder to be invariant to the specific
ordering of sequences in the input.

##### The input embedding

is formed by adding a learned embedding for each amino acid to the random
Fourier feature embedding of the column index as follows. First, the amino
acid token is embedded

$x_{i,j}^{aa}=W_{x_{i,j}},x_{i,j}^{aa}\in R^{d},W\in R^{K\times d}$

where $W$ is a matrix of learnable amino acid embeddings of dimension $d$, K
is the size of the vocabulary (22 in the case of 20 amino acids plus gap and
start/end tokens), and $x_{i,j}$ indicates the amino acid at position $i,j$ of
the MSA. Next, the column index is embedded

$x_{i,j}^{pos}=W^{pos}\cos(rj+b),$

$x_{i,j}^{pos}\in R^{d},r\in R^{d},b\in R^{d},W^{pos}\in R^{d\times d}$, where
$W^{pos}$ is a learnable matrix, $r$ is a random vector drawn from Normal(0,
1), and $b$ Is a random vector drawn from Uniform(0, 2$\pi$). The input
embedding is then formed by

$z_{i,j}=x_{i,j}^{aa}+x_{i,j}^{pos}$

##### The MSA encoder transformer layers

are composed of axial self-attention along the rows and columns of the MSA
matrix followed by a fully connected layer. Each operation is preceded by
layer normalization and uses residual connections (Algorithm 1).

Algorithm 1 Transformer MSA Encoder with Axial Attention

Input: MSA $z_{i,j}\in R^{d},i\in 1...N,j\in 1...M$

First, attention within each row

$z_{i,j}^{\prime}=LayerNorm(z_{i,j})$

$z_{i,j}^{\prime}=RowAttention(z_{i,j}^{\prime})$

$z_{i,j}=z_{i,j}+z_{i,j}^{\prime}$

Next, attention within each column

$z_{i,j}^{\prime}=LayerNorm(z_{i,j})$

$z_{i,j}^{\prime}=ColumnAttention(z_{i,j}^{\prime})$

$z_{i,j}=z_{i,j}+z_{i,j}^{\prime}$

Finally, the feed forward layer

$z_{i,j}^{\prime}=LayerNorm(z_{i,j})$

$h_{i,j}=GeLU(Linear(z_{i,j}^{\prime})),h_{i,j}\in R^{4d}$

$z_{i,j}^{\prime}=Linear(h_{i,j})$

$z_{i,j}=z_{i,j}+z_{i,j}^{\prime}$

Return: $z_{i,j}\in R^{d}$

Within row and within column attention are computed efficiently by calculating
the multi-headed attention operation over the rows or columns batch-wise. That
is, given a batch of intermediate MSA representations, Z, with dimensions
BxNxMxd, where B is the number of MSAs in the batch, we calculate per-row
self-attention by treating rows as part of the batch dimension, BNxMxd, and we
calculate per-column self-attention by treating columns as part of the batch
dimension, BMxNxd.

###  3.2 The Transformer MSA Decoder

##### The sequence decoder transformer layers

are composed of causal self-attention, cross-attention to the MSA
representations, and fully connected layers with layer normalization and
residual connection for each block (Algorithm 2). As is typical in transformer
decoders, the causal self-attention mask ensures that the output
representation for each position of the sequence in the decoder is only a
function of previous positions in the sequence and the MSA representations.
Before being passed into the decoder transformer layers, the target sequence
is embedded following the same scheme as the input embedding for the MSA along
the column dimension, $a_{k}$. When decoding, the target sequence is padded to
begin with a start token and end with a stop token.

Algorithm 2 Transformer Sequence Decoder with MSA Cross Attention

Input: MSA representations $z_{i,j}\in R^{d},i\in 1...N,j\in 1...M$, and
target sequence $a_{k}\in R^{d},k\in 1...L$

First, casual self-attention within the target sequence

$a^{\prime}_{k}=LayerNorm(a_{k})$

$a^{\prime}_{k}=CausalSelfAttention(a^{\prime}_{k})$

$a_{k}=a_{k}+a^{\prime}_{k}$

Next, cross attention against the MSA representations

$z_{i,j}^{\prime}=LayerNorm(z_{i,j})$

$a^{\prime}_{k}=LayerNorm(a_{k})$

$a^{\prime}_{k}=MSACrossAttention(a^{\prime}_{k},z_{i,j}^{\prime})$

Finally, the feed forward layer

$a^{\prime}_{k}=LayerNorm(a_{k})$

$h_{k}=GeLU(Linear(a^{\prime}_{k})),h_{k}\in R^{4d}$

$a^{\prime}_{k}=Linear(h_{k})$

$a_{k}=a_{k}+a^{\prime}_{k}$

Return: $a_{k}\in R^{d}$

In the MSA cross attention module, each position of the target sequence
attends over the complete MSA representations for each attention head (i.e.,
LxNxMxH, where H is the number of heads). We compute this efficiently by
flattening the MSA representations along the row and column dimensions such
that each MSA representation, $z_{i,j}$, is a single key in the cross
attention module.

Decoding the target sequence is performed using the decoder representations by
learning a transformation from $a_{k}$ to the probability distribution over
the $k+1$th token. We formulate this as a linear transformation of $a_{k}$
into a vector of dimension equal to the number of tokens followed by softmax
to give the probability of each token,
$p(y_{k+1}|a_{k})=p(y_{k+1}|y_{1...k},X)$.

###  3.3 Loss, hyperparameters, and training

##### Loss function.

We fit the parameters of the model to minimize the negative log-likelihood of
the target sequences conditioned on their family MSAs. Given $K$ MSA, target
sequence pairs, this loss is

$\mathcal{L}=-\sum_{i=1}^{K}\sum_{k=1}^{L^{i}}\log
p(y_{k}^{i}|y_{1...k-1}^{i},X^{i})$

##### Model hyperparameters.

Here, we train models with the following specific hyperparameters. We use 6
encoder and decoder layers each with hidden dimension ($d$) of 768. The MSA
encoder uses 12 attention heads and the decoder uses 8 attention heads.

##### Training details.

We train the model using the ADAM variant of stochastic gradient descent using
a linear ramp up, square root decay learning rate scheduler using a learning
rate of 0.0001 with 4,000 warmup steps. We use a total minibatch size of 256
spread over 16 GPUs using distributed training. Each GPU process minibatches
of size 1 with gradient accumulation over 16 steps to give the total effective
minibatch size. In order to reduce GPU RAM consumption, sequences and MSAs are
randomly sampled to a maximum length of 402 tokens during training.
Furthermore, MSAs are randomly downsampled to contain between 1 and 50
sequences. During training, we monitor the loss on a validation set of MSAs
and stop training when the validation loss stops decreasing.

##### Datasets.

We train and evaluate the model on full Pfam family alignments. In order to
evaluate the performance of the model on unseen families, we split the Pfam
sequences at the family level into 10,593 training, 563 validation, and 2,654
test families.

##  4 Results

###  4.1 Few shot protein generation of unseen protein families

We first examine the ability of our model to generate new proteins by
generalizing from few protein sequences representing a protein family to
unseen members of the family. We calculate the perplexity of sequences in the
validation set families given an increasing number of randomly selected
observed members, not including the target sequence. This allows us to
understand the ability of our model to extrapolate evolutionary landscapes
from a small number of observations. We compare our model with two widely used
protein statistical sequence models: position-specific scoring matrices
(PSSMs) and profile HMMs (pHMMs).

Table 1: Perplexities of generative models on the Pfam validation families as
the number of observed sequences increases. For each family, a target sequence
is sampled and then an MSA with the specified number of sequences is sampled
from the remaining sequences in the family. A PSSM or pHMM is then fit to the
MSA and used to calculate the perplexity of the target sequence or for our
model. For our model, the feed the MSA into the trained encoder and calculate
the perplexity of the target sequence using the trained decoder. We repeat
this 1,000 times for each family and report the average perplexity over trials
and families. Our model achieves lower perplexity at every MSA size and even
outperforms the profile HMM with 10x as many sequences.
| Number of sequences in MSA
---|---
Method | 1 | 10 | 25 | 50 | 100
PSSM | 14.30 | 8.43 | 7.17 | 6.64 | 6.20
pHMM | 12.30 | 6.65 | 5.85 | 5.51 | 5.32
Ours | 9.48 | 4.44 | 3.62 | 3.20 | 2.92

As expected, we observe that all models’ ability to generalize to unseen
family members improves as the number of observed family members increases
(Table 1). Furthermore, our model dramatically outperforms PSSMs and profile
HMMs across all MSA sizes. Remarkably, our model scales much better with
additional data, even outperforming PSSMs and pHMMs with 10x fewer sequences.
This demonstrates that our model learns how to extrapolate from small number
of sequences by better capturing evolutionary priors over sequences .

A major challenge in protein engineering is navigating the enormous search
space of possible sequence variants, because the space of sequence variants
increases exponentially with the number of sites . For example, if we consider
all 20 amino acids at 10 positions, the number of unique sequences is 2010
which is greater than ten trillion. At 65 sites, the space of possible
sequences exceeds the number of atoms in the universe. However, the vast
majority of these variants are not functional ($<$1% in typical mutagenesis
experiments). Therefore, homing in on only the space of viable protein
variants is critical for efficiently and feasibly searching sequence space.
Perplexity represents the number of amino acids that would need to be guessed
from uniformly to find the correct amino acid, therefore, it is the size of
the reduced alphabet learned by the model. On this basis, our model can
produce an enormous reduction in library size for protein engineering. Using
the 10 sites example, the pHMM perplexity of 5.3 yields a library size of 18
million which our model reduces to only 42,000, more than an order of
magnitude reduction in library size over the pHMM and about 8 orders of
magnitude better than random search. This improvement is even more extreme
when considering more sites for mutation.

Compared to the masked language models described, MSA2Prot offers exact
sampling through the use of a decoder. This sidesteps computationally
intensive Gibbs sampling.

###  4.2 MSA2Prot Log-Likelihood Predicts Variant Function

####  4.2.1 Point Mutants

We evaluate our model on the protein mutation datasets developed by [6]. These
datasets consist of mostly single mutant deep mutational scans, and a few
double mutant deep mutational scans.

As a baseline, we compare to an unconditional language model trained on the
same dataset and MSA2Prot predictions with only the wild-type sequence. We
also compare to to state-of-art methods for protein fitness prediction. These
include PSSMs, pHMMs, EVmutation, DeepSequence, ESM-1v, and MSATransformer.

PSSMs estimate the probability of each amino acid at a given position and
assume independence between different positions. EVmutation is a Potts model
that estimates the pairwise interaction terms between residues. pHMMs model
amino acids as being generated conditioned on a hidden state of the respective
MSA column. DeepSequence trains a generative model on a large multiple
sequence alignment for a given protein. ESM-1v and MSATransformer are masked
language models, where ESM-1v models protein sequences whereas MSATransformer
models multiple sequence alignments.

Table 2: Comparison of MSA2Prot to unsupervised protein fitness predictors.
New baselines include an unconditional language model and MSA2Prot with only
the wild-type sequence. Spearman Rho is averaged across 40 deep mutational
scanning datasets. HIS7 YEAST is omitted due to its large size.
Method | DMS Spearman Rho
---|---
MSA2Prot | 0.536
MSATransformer [21] | 0.548
ESM-1v (finetuned) [27] | 0.538
DeepSequence [6] | 0.513
ESM-1v (zero-shot) [27] | 0.508
EVMutation [28] | 0.506
PSSM | 0.460
MSA2Prot wildtype | 0.272
Unconditional LM | 0.091
Figure 2: Spearman Rho Across 40 DMS Datasets

Compared to the above methods, MSA2Prot displays state-of-the-art performance.
Comparisons with the unconditional language model as well as the wild-type
only MSA2Prot predictions indicate that the inclusion of the MSA is a key
driver of this performance. The averages across the datasets are shown in
Table 2, and a dataset-specific breakdown is shown in Figure 2 .

####  4.2.2 Higher Order Mutants

Given that 37 of above 40 DMS datasets consist of single mutants, we evaluate
MSA2Prot on a dataset of chorismate mutase sequences [29]. This dataset is
highly diverse, with sequence divergence in the test set ranging from 14
percent to 82 percent.

Table 3: Baseline spearman rho on chorismate mutate across methods. HMM, PSSM,
and Potts were fit to the MSA of natural sequences created by [29]. For
MSA2Prot, sequences ensembled with MSA size 100 and were reweighted with
sequences above 0.7 percent similarity considered to be neighbors. For MSA
Transformer, the best value across sequence reweighting with different
thresholds,random sampling, and the raw msa is reported.
Method | Spearman Rho
---|---
MSA2Prot | 0.41
MSATransformer | 0.12
HMM | 0.38
Potts | 0.41
PSSM | 0.39

We find that MSA2Prot outperforms MSATransformer substantially, as shown in
Table 3. This is likely because MSATransformer uses additive scoring to assess
the likelihood of higher order mutants, whereas MSA2Prot’s decoder is able to
exactly model epistasis.

Table 4: Spearman Rho for standard MSA, MSA filtered for variants above
fitness level 0.4, and MSA filtered for variants above fitness level 0.7. All
MSAs were reweighted with sequences above 0.7 percent similarity considered to
be neighbors. For the ensemble, MSAs of size 100 were sampled. MSA2Prot is
able to utilize negative samples, the 350 worst performing training sequences,
even when they have no predictive value on their own. MSA2Prot is also able to
generalize from the distribution of high-performing mutants.
| Fitness Level of MSA
---|---
Method | None | 0.4 | 0.7
Ensemble | 0.41 | 0.57 | 0.58
No Ensemble | 0.37 | 0.58 | 0.60
Ensemble w/ Neg Samples | 0.55 | 0.58 | 0.59
No Ensemble w/ Neg Samples | 0.50 | 0.59 | 0.58
PSSM | 0.39 | 0.50 | 0.50
HMM | 0.38 | 0.48 | 0.48
Negative Samples Only | 0.00 | 0.00 | 0.00

In addition, we find that MSA2Prot is able to generalize from a distribution
of high-performing mutants, as shown in Table 4. Spearman rho on the test set
dramatically improves after the training MSA is filtered to include only a few
hundred high-performing sequences. Lastly, we find that MSA2Prot is able to
harness low-performing variants to significantly improve accuracy. Although
these low-performing variants have no predictive power on their own, we use
them as negative samples by subtracting the likelihood of a sequence from the
low-performing MSA from the likelihood of a sequence from a standard MSA.

Often, combing the literature for a given protein will yield a list of high
and low performing mutants. However, given that experimental setups differ,
there may not be consistent fitness measurements. MSA2Prot is an ideal
candidate for this situation, given its ability to harness both high and low
performing variants without explicit functional measurements.

MSA2Prot also offers exact generation conditioned on multiple attributes.
Given a protein sequence, the probability distribution over the next residue
can be obtained by adding and re-normalizing the marginals of two MSAs, each
representing different attributes.

####  4.2.3 Indels

Table 5: Spearman Rho across methods on beta-lactamase indel dataset
Method | Spearman Rho
---|---
MSA2Prot | 0.46
DeepSequence | 0.45
HMM | 0.52

We evaluate MSA2Prot on [30]’s dataset of 262 deletions and 4422 insertions.
We benchmark against [6] and a HMM. Table 5 indicates that MSA2Prot achieves
comparable performance with [6], even though it does not require retraining.

###  4.3 Adaptive sampling with MSA2Prot

Table 6: Median log likelihood for test proteins with fitness greater than
0.4. Log likelihoods of test proteins were evaluated from an MSA with no
filtering, an MSA containing training sequences with fitness greater than 0.4,
and an MSA containing training sequences with fitness greater than 0.7. All
MSAs were reweighted with sequences above 0.7 percent similarity considered to
be neighbors. An input MSA with higher fitness leads to a higher likelihood on
high-performing test set variants
MSA Fitness | Log Likelihood
---|---
None | -124.36
0.4 | -112.52
0.7 | -108.92

In addition to offering high predictive accuracy through its log-likelihood,
MSA2Prot is able to generatively extrapolate from high-performing mutants. We
find that the high-performing sequences in the test set become several orders
of magnitude more likely as the training msa is filtered to include higher
performing variants. This result is shown by Table 6.

Figure 3: Change in Fitness on Chorismate Mutase Dataset Through Adaptive
Sampling. Sampled sequences replace MSA if they have a higher predicted
fitness, which is approximated by a Random Forest Regressor.Probability
distributions across each group of 500 sequences are shown

We also test MSA2Prot’s ability to adaptively sample high fitness variants
given black-box oracle,approximated by a Random Forest Regressor. We randomly
sample 100 sequences from the training MSA to form an initial MSA. We sample
sequences, and update the MSA if the regressor predicts the sampled sequence
has a higher fitness than the minimum fitness sequence in the MSA. The
results, shown in Figure 3, indicate that MSA2Prot is able to effectively
generate strong mutants. MSA2Prot initially generated sequences with a fitness
of 0.3, which is the 56th percentile of the training distribution. After
several thousand updates, MSA2Prot generated sequences with a fitness of 0.9,
which is in the 89th percentile of the training distribution.For comparison,
Gibbs sampling initially generated a sequence with a fitness of 0.03, which is
in the 39th percentile of the training distribution. After the same number of
function evaluations, Gibbs sampling generated a sequence with fitness 0.3,
which is in the 56th percentile of the training distribution. Thus, MSA2Prot
displayed considerably stronger performance. Compared to standard adaptive
sampling methods, MSA2Prot offers the benefit of not requiring training.
Instead, MSA2Prot relies on its few-shot generalization abilities. Thus,
MSA2Prot offers a computationally efficient alternative for protein sequence
design.

##  5 Conclusion

We present MSA2Prot, a transformer-based model for protein sequence generation
directly conditioned on related protein sequences represented as a multiple
sequence alignment. This model allows for efficient sampling and exact log-
likelihood evaluation for sequences including insertions and deletions as well
as substitutions without any need for training or fine-tuning on an MSA-by-MSA
basis. The MSA2Prot model generalizes well to unseen families and generates to
unseen family members with higher probability than other family generative
models, especially when MSAs are small. Furthermore, we demonstrate that this
model achieves state-of-the-art performance in variant function prediction on
diverse mutagenesis datasets. Finally, MSA2Prot unlocks new capabilities in
adaptive sequence generation by allowing the sampler to be focused only on
functional variants without expensive retraining.

## References

  * [1] Jiří Zahradník and Gideon Schreiber.  Protein engineering in the design of protein–protein interactions: Sars-cov-2 inhibitors as a test case.  volume 60, pages 3429–3435, 2021.  PMID: 34196543.
  * [2] Robert H Wilson, Elena Martin-Avila, Carly Conlan, and Spencer M Whitney.  An improved escherichia coli screen for rubisco identifies a protein-protein interface that can enhance CO2-fixation kinetics.  J. Biol. Chem., 293(1):18–27, January 2018.
  * [3] V J LiCata and G K Ackers.  Long-range, small magnitude nonadditivity of mutational effects in proteins.  Biochemistry, 34(10):3133–3139, March 1995.
  * [4] Tyler N Starr and Joseph W Thornton.  Epistasis in protein evolution.  Protein Sci., 25(7):1204–1218, July 2016.
  * [5] Kevin K Yang, Zachary Wu, and Frances H Arnold.  Machine-learning-guided directed evolution for protein engineering.  Nat. Methods, 16(8):687–694, August 2019.
  * [6] Adam J Riesselman, John B Ingraham, and Debora S Marks.  Deep generative models of genetic variation capture the effects of mutations.  Nat. Methods, 15(10):816–822, October 2018.
  * [7] Faruck Morcos, Andrea Pagnani, Bryan Lunt, Arianna Bertolino, Debora S Marks, Chris Sander, Riccardo Zecchina, José N Onuchic, Terence Hwa, and Martin Weigt.  Direct-coupling analysis of residue coevolution captures native contacts across many protein families.  Proc. Natl. Acad. Sci. U. S. A., 108(49):E1293–301, December 2011\.
  * [8] Sivaraman Balakrishnan, Hetunandan Kamisetty, Jaime G Carbonell, Su-In Lee, and Christopher James Langmead.  Learning generative models for protein fold families.  Proteins, 79(4):1061–1078, April 2011.
  * [9] Debora S Marks, Thomas A Hopf, and Chris Sander.  Protein structure prediction from sequence variation.  Nat. Biotechnol., 30(11):1072–1080, November 2012.
  * [10] Stefan Seemayer, Markus Gruber, and Johannes Söding.  CCMpred–fast and precise prediction of protein residue-residue contacts from correlated mutations.  Bioinformatics, 30(21):3128–3130, November 2014.
  * [11] Ronald M Levy, Allan Haldane, and William F Flynn.  Potts hamiltonian models of protein co-variation, free energy landscapes, and evolutionary fitness.  Curr. Opin. Struct. Biol., 43:55–62, April 2017.
  * [12] S F Altschul, T L Madden, A A Schäffer, J Zhang, Z Zhang, W Miller, and D J Lipman.  Gapped BLAST and PSI-BLAST: a new generation of protein database search programs.  Nucleic Acids Res., 25(17):3389–3402, September 1997.
  * [13] Michael Beckstette, Robert Homann, Robert Giegerich, and Stefan Kurtz.  Significant speedup of database searches with HMMs by search space reduction with PSSM family models.  Bioinformatics, 25(24):3251–3258, December 2009.
  * [14] Tristan Bepler and Bonnie Berger.  Learning protein sequence embeddings using information from structure.  In International Conference on Learning Representations, 2019.
  * [15] Tristan Bepler and Bonnie Berger.  Learning the protein language: Evolution, structure, and function.  Cell Syst., 12(6):654–669.e3, June 2021.
  * [16] Ethan C Alley, Grigory Khimulya, Surojit Biswas, Mohammed AlQuraishi, and George M Church.  Unified rational protein engineering with sequence-based deep representation learning.  Nat. Methods, 16(12):1315–1322, December 2019.
  * [17] Ahmed Elnaggar, Michael Heinzinger, Christian Dallago, Ghalia Rehawi, Yu Wang, Llion Jones, Tom Gibbs, Tamas Feher, Christoph Angerer, Martin Steinegger, Debsindhu Bhowmik, and Burkhard Rost.  ProtTrans: Towards cracking the language of lifes code through self-supervised deep learning and high performance computing.  IEEE Trans. Pattern Anal. Mach. Intell., PP, July 2021.
  * [18] Alexander Rives, Joshua Meier, Tom Sercu, Siddharth Goyal, Zeming Lin, Jason Liu, Demi Guo, Myle Ott, C Lawrence Zitnick, Jerry Ma, and Rob Fergus.  Biological structure and function emerge from scaling unsupervised learning to 250 million protein sequences.  Proc. Natl. Acad. Sci. U. S. A., 118(15):e2016239118, April 2021\.
  * [19] Jung-Eun Shin, Adam J Riesselman, Aaron W Kollasch, Conor McMahon, Elana Simon, Chris Sander, Aashish Manglik, Andrew C Kruse, and Debora S Marks.  Protein design and variant prediction using autoregressive generative models.  Nat. Commun., 12(1):2403, April 2021.
  * [20] Brian Hie, Ellen D Zhong, Bonnie Berger, and Bryan Bryson.  Learning the language of viral evolution and escape.  Science, 371(6526):284–288, January 2021.
  * [21] Joshua Meier, Roshan Rao, Robert Verkuil, Jason Liu, Tom Sercu, and Alexander Rives.  Language models enable zero-shot prediction of the effects of mutations on protein function.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
  * [22] Ali Madani, Bryan McCann, Nikhil Naik, Nitish Shirish Keskar, Namrata Anand, Raphael R. Eguchi, Po-Ssu Huang, and Richard Socher.  Progen: Language modeling for protein generation.  bioRxiv, 2020.
  * [23] John Ingraham, Vikas Garg, Regina Barzilay, and Tommi Jaakkola.  Generative models for graph-based protein design.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [24] David Brookes, Hahnbeom Park, and Jennifer Listgarten.  Conditioning by adaptive sampling for robust design.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 773–782. PMLR, 09–15 Jun 2019.
  * [25] Julienne M Mullaney, Ryan E Mills, W Stephen Pittard, and Scott E Devine.  Small insertions and deletions (INDELs) in human genomes.  Hum. Mol. Genet., 19(R2):R131–6, October 2010.
  * [26] David Brookes, Hahnbeom Park, and Jennifer Listgarten.  Design by adaptive sampling.  International conference on machine learning., 2019.
  * [27] Joshua Meier, Roshan Rao, Robert Verkuil, Jason Liu, Tom Sercu, and Alexander Rives.  Language models enable zero-shot prediction of the effects of mutations on protein function.  bioRxiv, 2021.
  * [28] Thomas A Hopf, John B Ingraham, Frank J Poelwijk, Charlotta P I Schärfe, Michael Springer, Chris Sander, and Debora S Marks.  Mutation effects predicted from sequence co-variation.  Nat. Biotechnol., 35(2):128–135, February 2017.
  * [29] William Russ, Matteo Figliuzzi, Christian Stocker, Pierre Barrat-Charlaix, and Peter Socolich.  An evolution-based model for designing chorismate mutase enzymes.  Science, 369(6502):440–445, July 2020.
  * [30] Courtney E. Gonzalez, Paul Roberts, and Marc Ostermeier.  Fitness effects of single amino acid insertions and deletions in tem-1 beta-lactamase.  Journal of Molecular Biology, 431(10):2320–2330, May 2019.

##  Appendix A DMS experiment

[6] compiled a set of 41 deep mutational scanning experiments, almost all
single-point mutations, and we used the MSAs provided by them. MSA filtering
is described section C below. HIS7 YEAST was omitted from the comparison due
to its large size of 496,000.

##  Appendix B Chorismate Mutase Dataset

[29] created a MSA of 1130 natural chorismate mutase enzymes. They sampled
1618 sequences from Potts models trained on this MSA. The MSA of natural
enzymes was used as a training MSA, and the 1618 variants were used as a test
set. Fitness was measured for sequences in the training and test set.

###  B.1 Spearman Rho Comparison Across Methods

For MSA2Prot,the MSA was ensembled by sampling groups of 100 sequences, where
a sequence’s probability of selection was inversely proportional to its number
of neighbors. Two sequences were considered neighbors if their sequence
similarity was greater than 0.7.The total number of groups of 100 was
determined by the number of effective sequences in the MSA. For
MSATransformer, the best spearman rho across sequence reweighting with
different thresholds,random sampling, and the raw msa is reported. Perplexity
was calculated using masked marginals and the wild-type sequence was chosen to
be E.coli. All other baselines were trained on the natural enzyme MSA without
any filtering. ESM-1v was omitted as a baseline due to the many insertions in
the MSA, relative to the wild-type.

###  B.2 Spearman Rho Comparison Across MSA Fitness

The training MSA was filtered to only include sequences above a threshold.
HMMs and PSSMs were trained on this MSA without any further filtering. If
ensembled, MSAs of size 100 were sampled from the training MSA, where a
sequence’s probability of selection was inversely proportional to its number
of neighbors. Two sequences were considered neighbors if their sequence
similarity was greater than 0.7. The total number of sequences sampled was
determined by the number of effective sequences in the MSA. For negative
sampling, an MSA was creating with the 350 worst-performing training
sequences. To generate spearman rho values, the predicted likelihood of a
sequence was the likelihood under the negative MSA subtracted from twice the
likelihood of the sequence under the standard MSA.

###  B.3 Log Likelihood Comparison Across MSA Fitness

Three MSAs were created, an MSA with no filtering, an MSA containing training
sequences with fitness greater than 0.4, and an MSA containing training
sequences with fitness greater than 0.7. Each of these MSAs were further
reweighted as described above, and the log likelihood of test proteins with
fitness greater than 0.4 was measured.

###  B.4 Adaptive Sampling

A random forest regression model with 1000 estimators was trained on the
training and test MSAs as a baseline. 100 sequences from the training MSA were
randomly sampled to form an MSA. Sequences from this MSA were sampled and
evaluated using the regressor. If a sampled sequence had a higher fitness than
the minimum fitness of the sequences in the MSA, then the minimum fitness
sequence in the MSA would be replaced by the sampled sequence.

The last 3 sequences from the adaptive sampling experiment are shown below.

LALGQLADERRDRDAQALALIDARQALAERVREVKRAAGAEVTDATRERAVLDRAEAGAVPQGEQNGRQGDSVLPELFAALYAQGRQAEALRARW

LHERLAIQVRMALDRELLFLIDARLTFAQQIGQIKRGKGAPVSDPSRQRVLQERAEANARAPQLAQPQLLVDLFERLFARSHTAERIQKRL

LELSYARQQQIDVDRRLLAMLAQRLPLAQQINNMKRSGGALVADPLRSRVIEERSNLMAATHSSDEMPFISDQFGRLFERVRNAERQARRR

##  Appendix C Beta-lactamase dataset

This deep mutational scanning dataset includes 4422 single insertions and 262
single deletions. An MSA was created using MMSeqs2, and was filtered and
ensembled as described above. The HMM was fit on the MSA, without any
filtering.

##  Appendix D Details of the MSA sampling strategies experiments

We experimented with several strategies for MSA selection. We varied MSAs by
size, ensembled MSAs, and combined MSAs with different sequence
similarities.We tested these strategies on a subset of five proteins from the
DeepSequence dataset. We report the average spearman rho across all of the
measurements for each protein.

We found that increasing MSA size, averaging predictions across multiple
sampled MSAs, and averaging predictions from MSAs sampled with different
sequence similarities all lead to improved performance. Thus, our final MSA
strategy consisted of averaging six MSAs of size 100. For each sequence
similarity filter, 1.0, 0.7, and 0.4, we sample a MSA twice.

For a given sequence similarity threshold, we removed sequences where the
percentage of amino acids that differ from the wild type protein was greater
than the threshold. A probability distribution over the remaining sequences
was constructed by weighting sequences in proportion to their hamming distance
from other sequences.

| Number of sequences in MSA |
---|---|---
Sequence Sim | 60 | 300 | 600 | 300 (ensem) | 600 (ensem)
1.0 | 0.33 | 0.35 | 0.36 | 0.36 | 0.43
0.7 | 0.31 | 0.34 | 0.37 | 0.40 | 0.38
0.4 | 0.33 | 0.38 | 0.39 | 0.33 | 0.38
All | 0.32 | 0.37 | 0.37 | 0.37 | 0.42


* (33) S. Cullen and M. Perelstein, “SN1987A constraints on large compact dimensions,” Phys. Rev. Lett. 83 (1999), 268-271 [arXiv:hep-ph/9903422 [hep-ph]]; V. D. Barger, T. Han, C. Kao and R. J. Zhang, “Astrophysical constraints on large extra dimensions,” Phys. Lett. B 461 (1999), 34-42 [arXiv:hep-ph/9905474 [hep-ph]]; C. Hanhart, D. R. Phillips, S. Reddy and M. J. Savage, “Extra dimensions, SN1987a, and nucleon-nucleon scattering data,” Nucl. Phys. B 595 (2001), 335-359 [arXiv:nucl-th/0007016 [nucl-th]]; C. Hanhart, J. A. Pons, D. R. Phillips and S. Reddy, “The Likelihood of GODs’ existence: Improving the SN1987a constraint on the size of large compact dimensions,” Phys. Lett. B 509 (2001), 1-9 [arXiv:astro-ph/0102063 [astro-ph]]; S. Hannestad and G. Raffelt, “New supernova limit on large extra dimensions,” Phys. Rev. Lett. 87 (2001), 051301 [arXiv:hep-ph/0103201 [hep-ph]]; S. Hannestad and G. G. Raffelt, “Stringent neutron star limits on large extra dimensions,” Phys. Rev. Lett. 88 (2002), 071301 [arXiv:hep-ph/0110067 [hep-ph]]; S. Hannestad and G. G. Raffelt, “Supernova and neutron star limits on large extra dimensions reexamined,” Phys. Rev. D 67 (2003), 125008 [erratum: Phys. Rev. D 69 (2004), 029901] [arXiv:hep-ph/0304029 [hep-ph]]; G. D. Starkman, D. Stojkovic and M. Trodden, “Large extra dimensions and cosmological problems,” Phys. Rev. D 63 (2001), 103511 [arXiv:hep-th/0012226 [hep-th]]; S. Hannestad, “Strong constraint on large extra dimensions from cosmology,” Phys. Rev. D 64 (2001), 023515 [arXiv:hep-ph/0102290 [hep-ph]].
  * (34) D. Atwood, C.P. Burgess, E. Filotas, F. Leblond, D. London and I. Maksymyk, “Supersymmetric large extra dimensions are small and/or numerous,” Phys. Rev. D 63 (2001), 025007 [arXiv:hep-ph/0007178 [hep-ph]]; C.P. Burgess, J. Matias and F. Quevedo, “MSLED: A Minimal supersymmetric large extra dimensions scenario,” Nucl. Phys. B 706 (2005), 71-99 [arXiv:hep-ph/0404135 [hep-ph]]; G. Azuelos, P. H. Beauchemin and C.P. Burgess, “Phenomenological constraints on extra dimensional scalars,” J. Phys. G 31 (2005) no.1, 1-19 [arXiv:hep-ph/0401125 [hep-ph]]; P. H. Beauchemin, G. Azuelos and C.P. Burgess, “Dimensionless coupling of bulk scalars at the LHC,” J. Phys. G 30 (2004), N17 [arXiv:hep-ph/0407196 [hep-ph]]; P. Callin and C.P. Burgess, “Deviations from Newton’s law in supersymmetric large extra dimensions,” Nucl. Phys. B 752 (2006), 60-79 [arXiv:hep-ph/0511216 [hep-ph]].
  * (35) L. Alvarez-Gaume and E. Witten, “Gravitational Anomalies,” Nucl. Phys. B 234 (1984), 269; M. B. Green and J. H. Schwarz, “Anomaly Cancellation in Supersymmetric D=10 Gauge Theory and Superstring Theory,” Phys. Lett. B 149 (1984), 117-122.
  * (36) S. Randjbar-Daemi, A. Salam, E. Sezgin and J. A. Strathdee, “An Anomaly Free Model in Six-Dimensions,” Phys. Lett. B 151 (1985), 351-356; M. B. Green, J. H. Schwarz and P. C. West, “Anomaly Free Chiral Theories in Six-Dimensions,” Nucl. Phys. B 254 (1985), 327-348; J. Erler, “Anomaly cancellation in six-dimensions,” J. Math. Phys. 35 (1994), 1819-1833 [arXiv:hep-th/9304104 [hep-th]]; J. H. Schwarz, “Anomaly - free supersymmetric models in six-dimensions,” Phys. Lett. B 371 (1996), 223-230 [arXiv:hep-th/9512053 [hep-th]].
  * (37) L. E. Ibanez and A. M. Uranga, “String theory and particle physics: An introduction to string phenomenology,” Cambridge University Press (2012).
  * (38) G.K. Leontaris and P. Shukla, “Stabilising all Kähler moduli in perturbative LVS,” [arXiv:2203.03362 [hep-th]].
  * (39) C.P. Burgess, J. M. Cline and M. Postma, “Axionic D3-D7 Inflation,” JHEP 03 (2009), 058 [arXiv:0811.1503 [hep-th]].
  * (40) A. Salam and E. Sezgin, “Chiral Compactification on Minkowski x S**2 of N=2 Einstein-Maxwell Supergravity in Six-Dimensions,” Phys. Lett. B 147 (1984), 47.
  * (41) H. Nishino and E. Sezgin, “Matter and Gauge Couplings of N=2 Supergravity in Six-Dimensions,” Phys. Lett. B 144 (1984), 187-192; H. Nishino and E. Sezgin, “The Complete $N=2$, $d=6$ Supergravity With Matter and Yang-Mills Couplings,” Nucl. Phys. B 278 (1986), 353-379.
  * (42) P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. (2020) 083C01.
  * (43) C. O’Hare, “cajohare/axionlimits: Axionlimits,” (2020), URL https://doi.org/10.5281/ zenodo.3932430.
  * (44) C.P. Burgess, J. J. H. Enns, P. Hayman and S. P. Patil, “Goldilocks Models of Higher-Dimensional Inflation (including modulus stabilization),” JCAP 08 (2016), 045 [arXiv:1605.03297 [gr-qc]].
  * (45) G. Efstathiou and J. R. Bond, “Isocurvature cold dark matter fluctuations,” Mon. Not. Roy. Astron. Soc. 218 (1986) no.1, 103-121; H. Kodama and M. Sasaki, “Evolution of Isocurvature Perturbations. 1. Photon - Baryon Universe,” Int. J. Mod. Phys. A 1 (1986), 265; A. D. Linde and V. F. Mukhanov, “Nongaussian isocurvature perturbations from inflation,” Phys. Rev. D 56 (1997), R535-R539 [arXiv:astro-ph/9610219 [astro-ph]]; P. J. E. Peebles, “An isocurvature cdm cosmogony. I. a worked example of evolution through inflation,” Astrophys. J. 510 (1999), 523 [arXiv:astro-ph/9805194 [astro-ph]]; K. Enqvist, H. Kurki-Suonio and J. Valiviita, “Limits on isocurvature fluctuations from boomerang and MAXIMA,” Phys. Rev. D 62 (2000), 103003 [arXiv:astro-ph/0006429 [astro-ph]]; N. Bartolo, S. Matarrese and A. Riotto, “Adiabatic and isocurvature perturbations from inflation: Power spectra and consistency relations,” Phys. Rev. D 64 (2001), 123504 [arXiv:astro-ph/0107502 [astro-ph]]; P. Crotty, J. Garcia-Bellido, J. Lesgourgues and A. Riazuelo, “Bounds on isocurvature perturbations from CMB and LSS data,” Phys. Rev. Lett. 91 (2003), 171301 [arXiv:astro-ph/0306286 [astro-ph]]; F. Di Marco, F. Finelli and R. Brandenberger, “Adiabatic and isocurvature perturbations for multifield generalized Einstein models,” Phys. Rev. D 67 (2003), 063512 [arXiv:astro-ph/0211276 [astro-ph]]; C. Gordon and W. Hu, “A Low CMB quadrupole from dark energy isocurvature perturbations,” Phys. Rev. D 70 (2004), 083003 [arXiv:astro-ph/0406496 [astro-ph]]; C. T. Byrnes and D. Wands, “Curvature and isocurvature perturbations from two-field inflation in a slow-roll expansion,” Phys. Rev. D 74 (2006), 043529 [arXiv:astro-ph/0605679 [astro-ph]]; Z. Lalak, D. Langlois, S. Pokorski and K. Turzynski, “Curvature and isocurvature perturbations in two-field inflation,” JCAP 07 (2007), 014 [arXiv:0704.0212 [hep-th]]; M. Beltran, J. Garcia-Bellido and J. Lesgourgues, “Isocurvature bounds on axions revisited,” Phys. Rev. D 75 (2007), 103507 [arXiv:hep-ph/0606107 [hep-ph]].
  * (46) Y. Akrami et al. [Planck], “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020), A10 [arXiv:1807.06211 [astro-ph.CO]].
  * (47) T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, “Modified Gravity and Cosmology,” Phys. Rept. 513 (2012), 1-189 [arXiv:1106.2476 [astro-ph.CO]].
  * (48) Y. Aghababaie, C.P. Burgess, S. L. Parameswaran and F. Quevedo, “SUSY breaking and moduli stabilization from fluxes in gauged 6-D supergravity,” JHEP 03 (2003), 032 [arXiv:hep-th/0212091 [hep-th]].
  * (49) C.P. Burgess and L. van Nierop, “Technically Natural Cosmological Constant From Supersymmetric 6D Brane Backreaction,” Phys. Dark Univ. 2 (2013), 1-16 [arXiv:1108.0345 [hep-th]]. C.P. Burgess, L. van Nierop and M. Williams, “Gravitational Forces on a Codimension-2 Brane,” JHEP 04 (2014), 032 [arXiv:1401.0511 [hep-th]]. C.P. Burgess, R. Diener and M. Williams, “The Gravity of Dark Vortices: Effective Field Theory for Branes and Strings Carrying Localized Flux,” JHEP 11 (2015), 049 [arXiv:1506.08095 [hep-th]]. C.P. Burgess, R. Diener and M. Williams, “EFT for Vortices with Dilaton-dependent Localized Flux,” JHEP 11 (2015), 054 [arXiv:1508.00856 [hep-th]]. C.P. Burgess, R. Diener and M. Williams, “A problem with $\delta$-functions: stress-energy constraints on bulk-brane matching (with comments on arXiv:1508.01124),” JHEP 01 (2016), 017 [arXiv:1509.04201 [hep-th]]. C.P. Burgess, R. Diener and M. Williams, “Self-Tuning at Large (Distances): 4D Description of Runaway Dilaton Capture,” JHEP 10 (2015), 177 [arXiv:1509.04209 [hep-th]].
  * (50) C.P. Burgess, R. C. Myers and F. Quevedo, “On spherically symmetric string solutions in four-dimensions,” Nucl. Phys. B 442 (1995), 75-96 [arXiv:hep-th/9410142 [hep-th]].
  * (51) K. S. Choi and J. H. Park, “Post Newtonian Feasibility of the Closed String Massless Sector,” [arXiv:2202.07413 [hep-th]].
  * (52) J. E. Moody and F. Wilczek, “New Macroscopic Forces?,” Phys. Rev. D 30 (1984), 130
  * (53) H. Georgi and L. Randall, “Flavor Conserving CP Violation in Invisible Axion Models,” Nucl. Phys. B 276 (1986), 241-252; I. B. Khriplovich, “Quark Electric Dipole Moment and Induced $\theta$ Term in the Kobayashi-Maskawa Model,” Phys. Lett. B 173 (1986), 193-196.
  * (54) R. J. Crewther, P. Di Vecchia, G. Veneziano and E. Witten, “Chiral Estimate of the Electric Dipole Moment of the Neutron in Quantum Chromodynamics,” Phys. Lett. B 88 (1979), 123 [erratum: Phys. Lett. B 91 (1980), 487]
  * (55) R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977), 1440-1443; S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978), 223-226 F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978), 279-282
  * (56) P. Di Vecchia and G. Veneziano, “Chiral Dynamics in the Large n Limit,” Nucl. Phys. B 171 (1980), 253-272; G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, “The QCD axion, precisely,” JHEP 01 (2016), 034 [arXiv:1511.02867 [hep-ph]].
  * (57) C. P. Burgess, M. W. Horbatsch and S. P. Patil, “Inflating in a Trough: Single-Field Effective Theory from Multiple-Field Curved Valleys,” JHEP 01 (2013), 133 [arXiv:1209.5701 [hep-th]].
  * (58) C. P. Burgess, H. M. Lee and M. Trott, “Power-counting and the Validity of the Classical Approximation During Inflation,” JHEP 09 (2009), 103 [arXiv:0902.4465 [hep-ph]]; P. Adshead, C. P. Burgess, R. Holman and S. Shandera, “Power-counting during single-field slow-roll inflation,” JCAP 02 (2018), 016 [arXiv:1708.07443 [hep-th]].
  * (59) T. D. Cohen, R. J. Furnstahl, and D. K. Griegel, “Quark and gluon condensates in nuclear matter,” Phys. Rev. C45 (1992) 1881–1893.
  * (60) A. Arvanitaki and S. Dubovsky, “Exploring the String Axiverse with Precision Black Hole Physics,” Phys. Rev. D 83 (2011), 044026 [arXiv:1004.3558 [hep-th]]; A. Arvanitaki, M. Baryakhtar and X. Huang, “Discovering the QCD Axion with Black Holes and Gravitational Waves,” Phys. Rev. D 91 (2015) no.8, 084011 [arXiv:1411.2263 [hep-ph]]; A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky and R. Lasenby, “Black Hole Mergers and the QCD Axion at Advanced LIGO,” Phys. Rev. D 95 (2017) no.4, 043001 [arXiv:1604.03958 [hep-ph]]; M. J. Stott, “Ultralight Bosonic Field Mass Bounds from Astrophysical Black Hole Spin,” [arXiv:2009.07206 [hep-ph]]; V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. Mcallister and M. J. Stott, “Superradiance Exclusions in the Landscape of Type IIB String Theory,” [arXiv:2011.08693 [hep-th]]; M. Baryakhtar, M. Galanis, R. Lasenby and O. Simon, “Black hole superradiance of self-interacting scalar fields,” Phys. Rev. D 103 (2021) no.9, 095019 [arXiv:2011.11646 [hep-ph]]; V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. McAllister and M. J. Stott, “Superradiance in string theory,” JCAP 07 (2021), 033 [arXiv:2103.06812 [hep-th]].
  * (61) A. Hook and J. Huang, “Probing axions with neutron star inspirals and other stellar processes,” JHEP 06 (2018), 036 [arXiv:1708.08464 [hep-ph]].
  * (62) J. Zhang, Z. Lyu, J. Huang, M. C. Johnson, L. Sagunski, M. Sakellariadou and H. Yang, “First Constraints on Nuclear Coupling of Axionlike Particles from the Binary Neutron Star Gravitational Wave Event GW170817,” Phys. Rev. Lett. 127 (2021) no.16, 161101 [arXiv:2105.13963 [hep-ph]].
  * (63) B. P. Abbott et al. [LIGO Scientific and Virgo], “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119 (2017) no.16, 161101 [arXiv:1710.05832 [gr-qc]].
  * (64) J. Khoury and A. Weltman, “Chameleon fields: Awaiting surprises for tests of gravity in space,” Phys. Rev. Lett. 93 (2004), 171104 [arXiv:astro-ph/0309300 [astro-ph]].
  * (65) C. Burrage and J. Sakstein, “Tests of Chameleon Gravity,” Living Rev. Rel. 21 (2018) no.1, 1 [arXiv:1709.09071 [astro-ph.CO]].
  * (66) P. Brax et al.. To appear.