Year: 2021

Abstract:
The design of novel protein sequences is providing paths towards the development of novel therapeutics and materials. At the forefront is the challenging field of de novo protein design, which looks to design protein sequences unlike those found in nature using general design methodologies. In this work, we develop a tool for de novo design, based on a deep generative sequence model, that rapidly samples novel protein sequences with diverse and ordered structures. To build this tool we develop a framework, called DARK, that trains the underlying generative model on an iteratively expanding set of synthetic sequences. The resulting model generalizes where models trained on natural sequences struggle and greatly improves on the efficiency of comparable sampling-based approaches. We further show how it can generate high quality candidates for de novo design problems and aid in the development of further novel design methods, in all, providing another step, amongst others, towards truly automated and intelligent protein design.

heading: 1 INTRODUCTION
text: Generative modelling is beginning to be used for the task of designing protein molecules, a problem which offers potential solutions to a vast number of medical (Chevalier et al., 2017; Silva et al., 2019) and industrial challenges (King et al., 2012; Wang et al., 2021). Computational protein design methods look to efficiently generate large numbers of candidate sequences, prior to laboratory validation, that are confidently predicted to have stable and ordered structures together with prespecified structural and functional attributes. Recent generative modelling studies have looked to improve on the high time and computational cost of contemporary design methods, which rely on Monte Carlo sampling and energy function-based physics simulations (Huang et al., 2016).
Here, our aim is to develop a tool, based on a deep generative model, that can rapidly generate sequences with stable and ordered structures. Furthermore, we build this tool for the important task of de novo protein design. Computational de novo protein design is arguably the most promising and general approach to protein design, but also the most difficult (Huang et al., 2016; Korendovych & DeGrado, 2020). As an approach, it looks to design proteins unlike those seen in nature, with desired structural or functional attributes. Crucially, this is done without using any information from a preexisting protein as a starting point, scaffold, or guide.1 Instead, the designer must rely on general design methods with little, if any, provided information to design suitable sequences with suitable predicted structures (Woolfson, 2021; Marcos et al., 2018; Vorobieva et al., 2021). This presents a remarkably challenging design setting, evident in that de novo designs have yet to be successfully generated for a number of common but biotechnologically important folds, like immunoglobins.
A tool that quickly generates sequences, unlike those seen in nature, with a diverse range of ordered structures is of great value as it can be used to both rapidly generate candidates directly for design tasks as well as provide a basis for developing further novel design methods. In this work, we develop deep generative sequence models that we show as satisfying this specification. Unlike previous generative protein design approaches that typically rely on natural sequences and structurebased conditioning information, we build an approach for de novo design by developing a framework for training unconditional autoregressive language models on synthetic sequences. We call
1The opposite of such an approach is protein engineering, an important but distinctly different field we do not cover here. However, for clarity, an example of protein engineering in a machine learning context would be fine-tuning a pre-trained language model on existing proteins with the desired attributes.
this framework Design in Areas of Restricted Knowledge (DARK) 2. DARK iteratively generates large datasets of synthetic sequences that are optimized in an unsupervised way to have unspecified ordered structures that are highly likely under a supervised protein structure predictor (Figure 1). We build this tool as an unconditional model to avoid limiting its potential applications and also to show that it is possible to learn general relationships between sequence and structure without the aid of conditioning information.
We evaluate our approach by first demonstrating that DARK models satisfy our criteria for a de novo design tool. In particular, we show that the final model, DARK3, generates novel sequences with ordered structures as judged by state-of-the-art structure predictor AlphaFold (V2; Jumper et al. (2021)). We also test the structural generalization of our approach using a stringent unseen structure-based test set and find DARK models performs well. We demonstrate the applications of DARK3 as our tool with an example of designing a sequence with fold commonly used as a de novo scaffold for grafting functional sites. Finally, we show that DARK models enable the development of new methods by developing AlphaFold refinement, a novel and efficient approach for producing high confidence design candidates with AlphaFold. We also use it to de novo design a sequence with a high confidence immunoglobin fold.
In summary, we make the following contributions: (1) We show that unconditional generative models of protein sequences can learn distributions that capture general structure information by learning from synthetic sequences. (2) We propose a novel framework, DARK, for efficiently training deep generative models on synthetic protein sequences. (3) We provide the final DARK3 model as a flexible and fast tool for de novo design that generates novel sequences with diverse ordered structures. (4) We demonstrate AlphaFold refinement, a novel way to efficiently de novo design high confidence sequence candidates using AlphaFold and DARK.

heading: 2 RELATED WORK
text: Here, we discuss related work applying generative models to protein design as well as those that have focused on de novo design, and broader contemporary protein design methods.
Generative sequence modelling in protein design There have been a variety of machine learning approaches (Killoran et al., 2017; Wang et al., 2018; Norn et al., 2021) including generative modelling (Anand & Huang, 2018; Sabban & Markovsky, 2020; Linder et al., 2020) in protein design and protein engineering (Yang et al., 2018; Sinai & Kelsic, 2020). Unconditional and conditional generative models have been explored in a variety of problems adjacent to protein design, such as variant prediction (Riesselman et al., 2018; Shin et al., 2021). Other examples are recent work using
2This is inspired by the huge ‘dark’ areas of the sequence space that are unexplored by nature (Taylor et al., 2009; Perdigão et al., 2015).
language models to learn representations for downstream supervised tasks (Alley et al., 2019; Bepler & Berger, 2019; Elnaggar et al., 2020; Rao et al., 2021).
Most relevant to the work presented here is the set of models that have been developed to explore design using generative models conditioned on structural information. Although a range of architectures have been used, these models vary primarily in the representations of structural information used for conditioning. Amongst others, models have conditioned on low information descriptions of protein folds (weak conditioning) (Greener et al., 2018; Karimi et al., 2020), secondary structure labels (Singer et al., 2021), and course descriptions of protein backbones (Ingraham et al., 2019; Strokach et al., 2020; Cao et al., 2021).Some methods train proteins of one structural or functional type, effectively conditioning on that specific type (Yu & Buehler, 2020; Shin et al., 2021).
Conditioning out structure provides a way to address an inherent challenge to working with natural protein sequences. The vast majority of known natural proteins have either a weak signal for structure or it is mixed in with other information, such as functional constraints and evolutionary drift. Unconditional generative models, trained on natural sequences, fail to generalize to sequence with structures beyond those represented by their training sets because of this lack of structure signal (Ingraham et al., 2019). This is also well understood when viewed from the lens of protein structure prediction. Even highly accurate structure predictions models like AlphaFold (Jumper et al., 2021) rely on the information contained in aligned sets of sequences similar to the query sequence to achieve accurate predictions for natural sequences. With only individual natural sequences, performance tends to be poor (Moffat & Jones, 2021; Xu et al., 2021).
Generative Models and de novo protein design A generative sequence model for de novo design needs to be highly general, learning general sequence and structure relationships. All considered, we find it unsurprising that there has been no clear and unambiguous demonstration of a deep generative model, without any aid from conditioning, that generates sequences unlike natural sequences and, crucially, with confidently predicted ordered structures. Here, we demonstrate such a model. We achieve this by learning from syntheticsequences with a strong structure signal.
Previously de novo designed sequences are known to have extremely strong structure signals, but just over 100 de novo designed sequences (Woolfson, 2021) have lab-determined structures in the Protein Data Bank (PDB) (Burley et al., 2021). Thus, we find an alternative approach which is to generate synthetic sequences by leveraging a sequence optimization objective from the recent trDesign method (Anishchenko et al., 2020). Like contemporary approaches (Leaver-Fay et al., 2011), it is primarily a Monte Carlo sampling method, however in contrast, it slowly samples from a trained supervised structure predictor, using simulated annealing to maximize the structure likelihood of random sequences. The final sequences were shown to be extremely similar to pre-existing de novo designed sequences in their attributes and to have a strong structure signal.
Contemporary de novo protein design Although recent nonparametric methods have been developed (Zhou et al., 2020), by far the most historically successful de novo design approach is the Rosetta framework (Leaver-Fay et al., 2011), which is the recognized gold-standard (Huang et al., 2016). The Rosetta approach constitutes a complex computational pipeline that uses a variety of different protocols to iteratively sample viable structures and sequences using physics-based force-field simulations and Monte Carlo search procedures (Vorobieva et al., 2021; Marcos et al., 2018) making it lengthy and computationally expensive. New techniques that can rapidly produce quality samples, such as the work presented here, offer the potential to greatly reduce these costs. For further discussion regarding de novo design, we point the reader to the following reviews and retrospectives by Huang et al. (2016), Korendovych & DeGrado (2020), and (Woolfson, 2021).

heading: 3 LEARNING DARK MODELS

heading: 3.1 BACKGROUND
text: Let x ∈ V L be a protein sequence of length L over a vocabulary V of discrete amino acids (|V | = 20). Each amino acid is textually represented by a one-letter code. Each sequence x maps to a corresponding molecular structure y in three-dimensional space. Protein structure is described by the 3D coordinates of each atom, however more coarse-grained representations are commonly used. Here, we use the discrete representation of structure known as distograms (Senior et al., 2020). When using AlphaFold we use its predicted distograms and atomic coordinates, depending on the task (Jumper et al., 2021). Distograms represent a structure as a [L,L,D] tensor, comprising Ddimensional one-hot encodings of a discretized distance measure between residue pairs. We do not consider the amino acid Cysteine for experimental reasons when generating synthetic sequences3, however, for consistency we allow all trained models to predict it. For simplicity, we assume a fixed length of L = 100, however our approach easily generalizes to variable sequence lengths as it is based on autoregressively factorized probability models.
We are interested in the set of sequences that are unlike those in nature xd ∈ V L with structures that are folded in some ordered and regular state yd. We would like to sample from their joint distribution, which is convenient to break down, by chain rule, to p(xd, yd) = p(yd|xd)p(xd), a common way to consider the protein design problem. The first term p(yd|xd) is often thought of as protein structure prediction and so it is convenient to model it with some accurate proxy regression model pβ(yd|xd), being a trained supervised structure predictor, with parameters β ∈ B. We refer to this proxy regression model as an oracle. In this work, we focus on modeling the prior, p(xd) as the tool we desire, which, when combined with AlphaFold, provides a rapid means to produce high quality sequence and structure samples for de novo design problems.

heading: 3.2 AN UNSUPERVISED DESIGN LOSS FUNCTION
text: Recently, Anishchenko et al. (2020) demonstrated a method of sampling sequences for de novo design which, using our own notation, also looks to maximize p(yd|xd)p(xd), using the trRosetta structure predictor as an oracle (Yang et al., 2020). This samples sequences using simulated annealing of random sequences to optimize a combined proxy objective for both terms, ultimately generating sequences with the structure signal we are interested in. Here we describe the proxy objective in our own notation and then how we use it. First, p(yd|xd) is approximated using the Kullback-Leibler divergence (DKL) between an oracle and a trained background network p (yd) which approximates a prior over ordered structures.
arg max xd log pβ(yd|xd) ≈ arg max xd DKL
( pβ(yd|xd)‖p (yd) ) (1)
The background network is trained by taking a model with the same architecture and dataset as pβ(yd|xd) and training it with the sequence inputs replaced by Gaussian noise4 (Appendix A.1). The final objective also includes a term for approximating a weak sequence prior by taking the
3This is common in design studies, in part because it significantly complicates otherwise straightforward recombinant protein expression in bacteria, the most relevant factor for us (See Appendix A.7).
4In practice, for a given length, the background model distograms can be calculated once, taking an average over an arbitrary number of predictions (e.g. 100), and then saved to be used for all IG-score calculations.
negative DKL between the sequence’s normalized unigram frequencies p(xunid ) and the unigram frequencies of proteins with high resolution crystal structures in the PDB, p(xPDB). For simplicity, we refer to this objective as the Information Gain score (IG-score):
IG-score = arg max xd
( DKL ( pβ(yd|xd)‖p (yd) ) −DKL ( p(xunid )‖p(xPDB) )) (2)
Here, we implement our own version of trDesign, termed DARK0, and correspondingly our own version of the IG-score, using a version of the recent DMPfold2 model, chosen for its fast inference, for the oracle and background network (Kandathil et al., 2020), which is modified so that it produces distograms instead of atomic coordinates (Appendix A.1). DARK0 is slow, it takes an average of 11.4 minutes to generates a sample. Given the end point of protein design is in the lab, a performant design model needs to provide as many high confidence candidates as possible to justify the cost of validation as well as combat very high attrition rates. A method like trDesign is too slow, and has less potential future adaptions, to be used as the tool we desire. Instead we show how we can use IG-score to learn generative models that sample similar sequences orders of magnitude faster with many potential applications. To kickstart the approach described in the next section, we use DARK0 (40,000 simulated annealing steps) to sample a small number of sequences, 15K, that we call the seed samples. It is also used to sample two sets of 950 sequences to be used as a validation and test set for models in the next section (See Appendix B.1 for further details). After this DARK0 is discarded as an approach, and the only aspect of trDesign that we build on is our own version of the IG-score; the rest of what is described here is our own construction.

heading: 3.3 BUILDING THE DARK FRAMEWORK FOR DESIGN
text: Algorithm 1 Design in Areas of Restricted Knowledge (DARK) Input: Oracle model pβ(yd|xd) Input: Iterations N Input: Sample sizes at each iteration m0:N = {m0, ...,mN}
1: xs0 ← Generate m0 seed samples from pβ(yd|xd) 2: for n = 1, 2, ..., N do 3: pγn(xs)← pγn(xs|xsn−1), Learn a new sequence model 4: xs∗ ← Ancestral sampling of pγn(xs) . mn new samples 5: xs∗ ← Rapid optimization of xs∗ IG-score using the oracle 6: xsn ← [xs∗ , xsn−1 ] Combine optimized samples with existing samples 7: end for
DARK (Algorithm 1) tackles the primary limitation on effectively learning models, which is having suitable sequences to train and test against. The first step is to take some initial set of samples xs0 to train an initial sampling model pγ1(xs0) with parameters γ1 ∈ Γ. Here we use the aforementioned seed samples, but this could be random sequences or sequences from DARK0-Grad, an approach we developed that backpropagates the negative IG-score to generate samples (Appendix A.2), in a fashion similar to (Killoran et al., 2017). We use the seed samples out of convenience as they have the highest score, and we reasoned that they would require less iterations by comparison. The sampling models are all generative sequence models that can be sampled. After training pγ1(xs0), it is used to generate a large set of samples xs∗. Instead of then training a new sampling model on xs∗, each sequence in xs∗ is quickly optimized to improve the IG-score.
This is done by a quick greedy hill-climb on the sequence for a set number of steps (Line 5 in Algorithm 1), which we refer to as refinement. In each step, a random position is mutated. If it improves the IG-score then it is kept, otherwise it is discarded. We use 3000 steps, a fraction of the steps done in DARK0, as this was found to improve the average IG-score of xs∗ to approximately equal to the IG-score of xs0 . For simplicity, this is fixed for all steps. The result of optimizing xs∗ is then combined with xs0 to make xs1 . In the second iteration, xs1 is used to train a new sampling model pγ2(xs1). After some number of iterations N , the final sampling model is viewed as a strong prior over xd, also making it a powerful generative model that can be used for de novo design tasks. Here we perform 3 iterations with training set sizes of 15K, 100K, and 500K, where the difference is sampled and refined between iterations. These sample sizes are arbitrary and they were entirely dictated by computational resources available at the time. Three iterations were performed

heading: Language Model PPL Test ↓ IG-score ↑ IG-score (Ref) ↑ Sample Speed (s) ↓
text: as we believed it sufficient to demonstrate DARK’s effectiveness. We note that the iterative aspect of DARK shares similarities to methods like Estimation of Distribution Algorithms (EDAs), a common approach in model based optimization (MBO) (Bengoetxea et al., 2001; Larrañaga & Lozano, 2001; Brookes et al., 2020). Viewed through this lens, it suggests a number of ways that DARK can be adapted and potentially improved on, which we leave for future work (Kumar & Levine, 2020; Trabucco et al., 2021).
Deep generative architectures In DARK, we use deep autoregressive language models for learning protein sequence distributions. Specifically, we use a standard Transformer decoder architecture (Vaswani et al., 2017; Radford et al., 2019) which has been used extensively to train state-of-the-art language models capable of generating high-quality synthetic sequence samples. In the first iteration of DARK we use a small decoder model, termed DARK1, with 4 layers, 4 heads, and a feed-forward size of H = 128. To affirm the choice of self-attention based architecture, we compare it to a 1- layer LSTM with a hidden dimension H = 128. We also include results for a variant of DARK1 referred to as DARK1-Adversarial, that uses DMPfold2 as an adversarial regularizer during training (Appendix A.3). After the first iteration, the second and third iteration models, DARK2 and DARK3 respectively, use a decoder with 12 layers, 12 heads, and H = 768. The parameters are increased to coincide with the order of magnitude change in the size of the training set from DARK1 to DARK2. To provide a even comparison to models trained on natural sequences, we train a model, termed UF50, on 50M natural sequences from the UniRef50 sequence database (The UniProt Consortium, 2021) using the same architecture as DARK3 (Appendix A.5). For all trained models we perform early stopping and a small amount of hyper-parameter optimization on the validation set.
Quantitatively measuring sample quality We define high-quality samples as those that are diverse, have sequences unlike any natural sequence, and are predicted to have a stable and ordered structure. In de novo design, the most important test of a candidate sequence is that it is confidently predicted to have a stable and ordered atomistic structure (Huang et al., 2016). We provide a direct measure of both confidence and order jointly using AlphaFold’s pLDDT score, being its confidence metric, which ranges from 0 to 100 (Jumper et al., 2021). This is a measure of confidence, and an indirect measure of order. It was found that low scores (pLDDT< 50) strongly suggest disorder and vice versa (Tunyasuvunakool et al., 2021). Thus, we use the proportion of samples with a pLDDT> 70 as a measure of high quality for the predicted structures, which we refer to as Good+ pLDDT, indicating “good or better” quality according to criteria in Tunyasuvunakool et al. (2021). When discussing results we use ‘ordered’ to mean ‘confidently predicted and ordered’. Using AlphaFold also provides a check to ensure that the model has not learned to produce adversarial examples against the original DMPfold2 oracle. We also measure the diversity of sequences and predicted structures using the estimated number of clusters in samples of sequences, and samples of predicted structures (See Appendix B.7 & B.8).

heading: Language Model Good+ pLDDT (%) ↑ Training Data
text: Measuring the structural generality of DARK in a strict setting We take two key approaches to providing a stringent measure of how well DARK generalizes with regards to structure. Both are enabled by using AlphaFold to generate all-atom structure predictions for all 100,000 sequences contained in the DARK2 training set, which covers the 15,000 seed samples and the 85,000 samples from the first iteration of DARK. We then assign each sample’s structure to a topology. This process is the repeated for the 950 sequence validation set; we ignore the 950 sequence test set for what is described here. The first approach we take is to carefully construct a strict training and test set split based on structure, from the 100K samples (See Appendix B.2 for exhaustive details). This has been recently suggested for machine learning-based protein design studies (Ingraham et al., 2019) and is considered a gold-standard approach in protein structure prediction (Chothia & Lesk, 1986; Söding & Remmert, 2011). Our strict test set contains no sequences with overlapping topologies in either the training set or validation set. We also remove any sequences from the validation and training set that are detected to be similar to those in the test set with the MMseqs2 search tool (Steinegger & Söding, 2017). We train and evaluate a model on these sets, termed DARK2-STRICT. In the limit of working with synthetic, we believe this constitutes a stringent test of generality. The second approach we take is to compare the number of unique topologies represented by the 15,000 seed samples and those represented by the 85,000 sequences from the first iteration of DARK to show that DARK iterations generate sequences with unseen structures.
3.4 USING DARK3 FOR DE NOVO DESIGN
Rapid de novo design with DARK We provide an example of how DARK can be used in a specific design task by setting the challenge of de novo designing a 4-helix bundle, which is commonly seen in design studies and in nature (Orengo et al., 1997). 4-helix bundles are also of interest as they are being used as de novo scaffolds for grafting in functional sites (Woolfson, 2021). To design a sequence, we use DARK and AlphaFold to generate a database of samples and predicted structures. We take the highest pLDDT 4-helix bundle from the dataset and then validate its structure prediction by predicting its structure using the Rosetta Abinitio Relax protocol without any homology information (Leaver-Fay et al., 2011) (this took 433 CPU hours). We generate 10,000 structure predictions (decoys), and rank them by lowest Rosetta all-atom score. For this example, we use the set of 1000 unrefined sequences from DARK3 as if it were a structure database, but this can be expanded to an arbitrary size in practice. 1000 sequences with predicted structures is equivalent to generating 28 sequences with an approach like trDesign by time. Due to its speed, DARK is able to provide a large number of potential candidates sequences for different design tasks like this one in an off-the-shelf way.
An novel but simple approach to de novo designing with AlphaFold To demonstrate that DARK models can be used in method development, we use DARK3 to develop a simple design method we call AlphaFold refinement (Figure 2A & 2B). This provides a solution to the problem of having a candidate sequence with a structure of interest to a designer but it has low confidence. AlphaFold refinement functions by bootstrapping the initial predicted structure as a target for performing a hill-climb based optimization of the DARK3 sampled sequence. More specifically, we minimize the Cross-Entropy between the AlphaFold predicted distogram and the one-hot distogram of the bootstrapped target structure. To provide a comparison to AlphaFold refinement, we construct a trDesign-based approach that optimizes an AlphaFold IG-score (AIG-score) using simulated annealing. AIG-score also provides a measure of confidence. As an example of its practical use, we use AlphaFold refinement to de novo design a sequence with a high confidence immunoglobin fold (Figure 3) (See Appendix B.9 for exhaustive details).

heading: 4 RESULTS
text: Here, we assess DARK3 as a tool and provide demonstrations of its use.

heading: 4.1 EVALUATING DARK MODELS AS A TOOL FOR DE NOVO PROTEIN DESIGN
text: Comparing DARK to a simulated annealing approach We use DARK3 without refinement as our final design tool and we find that it produces similar samples by IG-score to DARK0 (seed samples), our version of trDesign, but does so four orders of magnitude faster, 52 ms compared to 11 mins per example (Table 1). We also find that the iterative approach to DARK effective for improving performance between iterations and significantly reduces the required resources. Rounding up to the nearest day, if we were to re-perform DARK from nothing to having a trained DARK3 it would take 12 days when parallelized across ten V100 GPUs. Of that time, model training constitutes just over 3 days and only requires 1 GPU. To generate 500,000 examples with DARK0 would take ∼80 days across ten V100 GPUs . Compared to methods like trDesign, the sample speed of DARK models allows for the generation of large quantities of candidate sequences for different design tasks. With an approach like DARK’s, that relies on generative neural networks, there are also clear ways to introduce controlled generation (Dathathri et al., 2019; Keskar et al., 2019).
DARK generates samples with ordered structures We find that DARK3 reliably generates samples with ordered predicted structures (Table 2). Without any refining or conditioning, 32.3% of DARK3 samples have Good+ pLDDT scores. Given the IG-score of unrefined DARK3 samples is slightly lower than that of the DARK0 seed samples (Table 1), it is unsurprising that similar difference is present with the Good+ pLDDT (Table 2). Compared to an approach like trDesign, for a small decrease in Good+ pLDDT, DARK3 is able to generate sequences orders of magnitude faster, also giving it much greater flexibility for integrating and developing future approaches.
DARK samples are diverse We find that samples from DARK models have diverse sequences and predicted structures (Table 3 & 4). Across all iterations, clustering results suggest that only a small number of samples have similar sequences and structures. That said, clustering is not a very robust measure as it doesn’t indicate an exact degree of diversity. We also generate 5000 sequences from DARK3 and predict their structures to gain a further qualitative insight into the diversity of the structures represented in the samples. There appears to be a preference for all-α and all-β structures, but α & β proteins are still common. We include examples in Figure 5. Looking forward, future investigation is required to properly quantify and examine the structure distribution of sequences from the final model. That said, we do find sequences of design interest in the samples that we have generated. For example, sequences with topologies that have yet to be successfully de novo designed, like immunoglobin folds.
DARK samples are unlike natural sequences We find that samples from the DARK3 model are distinct from sequences in nature. A sequence search, using the MMseqs2 search tool (Steinegger & Söding, 2017) with a high sensitivity parameter (-s 7), against UniRef100 (The UniProt Consortium, 2021) of the 1000 evaluated sequence samples returned no hits for all sequences except one, which had 5 weak hits. Repeating the sequence search with a different tool, HHblits (Remmert et al., 2012), no natural sequences were found that met the HHblits default threshold value, being an E-value below 1× 10−3 (See Appendix B.7 for further details). This strongly suggests that DARK3 generates sequences that are unlike those found in nature. We find the opposite results for the natural sequence model. Repeating the MMseqs2 search with the 1000 samples from the UF50 model results in 41 samples having a combined total of 1.76 million unique hits. Given the UF50 model is trained on natural sequences this is unsurprising. All considered, training on synthetic sequences appears to provide a simple means to generate sequences unlike those in nature.
Learning from natural sequences leads to poor performance The UF50 model provides a direct comparison to having trained on natural sequences instead of synthetic sequences, and is similar to protein engineering focused approaches like ProGEN (Madani et al., 2020) and UniRep (Alley et al., 2019). As expected, it performs poorly by both perplexity and IG-score (Table 1) and underperforms compared to DARK3 by Good+ pLDDT (Table 2). Especially given the test set comprises synthetic sequences, from a naive perspective we would expect UF50 to assign high likelihoods to the test set.
Not only is its training set two orders of magnitude larger than DARK3’s, the test set sequences have a clear structure signal, It is also hyperparameter optimized on the synthetic validation set. For the Good+ pLDDT, it is also advantaged as many of the sequences in the AlphaFold training set are in its own training set. Ultimately, this result is entirely expected (See Section 2).
DARK models generalize to sequences with unseen structures We find that DARK generalizes to sequences with unseen structures. DARK2-STRICT, tested on on the strict test set (Section 3.3 & Appendix B.2), performs well and achieves a perplexity of 6.50, similar to the other DARK models and similarly outperforming the baselines (Table 1). In assigning high likelihoods to the sequences in the strict test set, DARK2-STRICT provides direct evidence that DARK models are learning general relationships between sequences and structure. This is further supported by how unseen topologies are sampled between DARK iterations (Section 3.3). Within the 100K predicted structures of the DARK2, we find that the initial 15K seed samples contains 421 unique topologies and the 85K set from the first iteration of DARK contains 1011 unique topologies, with an intersection of 373 between the two. As such, there are 638 topologies sampled during the first DARK iteration that were not present in the 15K training set. This accounts for 1476 sequences, of which 225 sequences (15.2%) have a pLDDT≥ 70. This provides direct evidence that the DARK iterations are expanding the training sets to include sequences with unseen and ordered structures, suggesting a general approach.
4.2 USING DARK3 FOR DE NOVO DESIGN
An example of de novo design with DARK For the task of designing a 4-helix bundle (Section 3.4), we find that DARK3 successfully produces a candidate, effectively off-the-shelf. The predicted structure of the chosen sequence has a high pLDDT of 91.4, similar to the pLDDT of existing de novo proteins (Appendix B.4). We also find that its AlphaFold structure is in agreement with the best ranked predicted structure by Rosetta, having a low root-mean-square deviation of 3.11Å (by α-carbons) compared to AlphaFold’s prediction (Figure 6 & 8). In all, this provides a demonstration of how models like DARK3, paired with AlphaFold, can rapidly provide de novo candidates for tasks like generating scaffolds for grafting in functional sites.
Generating de novo designed sequences with AlphaFold refinement and DARK We find AlphaFold refinement can reliably generate high confidence (pLDDT 91.0) sequence candidates, while maintaining their structures,. This is slightly better than the confidence (pLDDT 88.1) with which AlphaFold predicts the structures of recently released de novo designed sequences in the PDB (Figure 2C & 2E). We also find that AlphaFold refinement is efficient, typically converging by 7500 steps which is ∼10 hours on a consumer-grade graphics card (RTX 2080 Ti). This is significantly more efficient than approaches like de novo design with Rosetta which can take thousands of CPU hours. In contrast, we found the simulated annealing approach produced what we consider an adversarial attack on AlphaFold (Figure 2E & 4). For the task of de novo designing a sequence with an immunoglobin fold, we find that DARK with AlphaFold refinement is also successful, producing a final predicted structure with a high pLDDT of 93.1 (Figure 3), demonstrating how DARK can be used to generate candidates for important but yet to be de novo designed folds.

heading: 5 CONCLUSION
text: In this work, we show that it is possible to build a tool, in the form of DARK3, that can rapidly generate sequences with diverse and ordered structures. Looking forward, building a variable length training set and smoothing the topology distribution between iterations are clear first targets as well as exploring related approaches in MBO. Controlled generation is an especially useful route to explore and there is a clear way to approach it: AlphaFold structure predictions themselves can be used as pseudo-labels for training future conditional versions of DARK entirely based on synthetic sequences. The AlphaFold refinement approach can also be easily extended to the sub-problem of fixed-backbone protein design. Looking forward, we hope that DARK, and many others, are continuing to provide steps towards intelligent and truly automatic protein design.

heading: ACKNOWLEDGMENTS
text: Left blank for anonymity during double-blind peer review.

heading: REPRODUCIBILITY STATEMENT
text: In order to guarantee reproduciblity of this work, we provide a detailed appendix with specific implementation and methodological details which is sufficient to reproduce the results in this work. Following release of the paper, we will release both a code repository containing the source code for the DARK models as well as the associated scoring functions and analysis scripts. The training data and model parameters will be made available from our servers due to GitHub file size limits.

heading: APPENDIX
text: A IMPLEMENTATION AND TRAINING DETAILS

heading: A.1 ORACLE MODEL AND BACKGROUND MODEL
text: We use the DMPfold2 model for the pβ(y|x) oracle (Kandathil et al., 2020). We received the trained model as well as the training and validation data described in Kandathil et al. (2020) via correspondence with the authors. We note that the code and pre-print for a more recent but non-probabilistic version appears to be available5, but we refer here to the CASP14 version of DMPfold2 described in the first pre-print (Kandathil et al., 2020), which outputs distograms.
DMPfold2 was chosen as an oracle as it reports favourable performance and, importantly, its ability to perform fast inference directly from an input alignment, or in our case, a single sequence. It uses a combination of two GRUs (Cho et al., 2014) followed by a residual network to predict distograms from a multiple sequence alignment. The input is of shape [A,L, 22] where the A dimension is the number of sequences in the alignment. First, a stack of unidirectional GRUs reduces A to a fixed size, resulting in an output array sized [L, 512], that is fed as input to a stack of bidirectional GRUs that produces the final representation of the input alignment. This representation is then fed to a convolutional ResNet, which produces the final outputs (see below).
We refactored the code for the DMPfold2 model in order to perform batched inference on individual sequences. We find that the IG-score for a single sequence (L = 100) can be calculated in approximately 0.06s. For a batch of 100 sequences this reduces to 0.01s per sequence.
Distogram outputs from DMPfold2 From a single input sequence of length L, DMPfold2 produces an output tensor of size [L,L, 104]. The 104D third axis is split into four separate discrete distributions which are individually softmaxed, being sized [L,L, 2] and three of [L,L, 34]. The output tensor contains predictions of residue-residue hydrogen bonding, residue-residue β-carbon (Cβ) distance distributions (distograms), and backbone φ and ψ torsion angle distributions. See Kandathil et al. (2020) for further details.
DMPfold2 Background Model To implement the IG-score we train a background network, p (yd), as per Anishchenko et al. (2020), using the DMPfold2 architecture and author-provided training data. This model is distinguished with the subscript . Training the background model constituted replacing the input sequences during training with 64D Gaussian noise, and altering the input layer accordingly, but otherwise following the training procedure described in the DMPfold2 paper. We pre-calculate the background values for a given sequence length by feeding noise to the network 100 times and taking the average of the predictions, leaving a final tensor of size [L,L, 104] containing the resulting predictions. Although this can be calculated on the fly, for simplicity we pre-calculate and save them once. As we focus on L = 100, this amounts to one saved array of size [100, 100, 104].
Using DMPfold2 to calculate the IG-score The IG-score is calculated by first predicting the output tensor for a sequence with DMPfold2. Next, DKL is calculated between each of the four distograms in the tensor and their corresponding pre-calculated background versions. This results in four [L,L, 1] tensors, one for each distogram, containing the scalar DKL values. For each tensor, the mean DKL is taken across both dimensions, resulting in four scalars. These four DKL values are added together with equal weighting. Finally, the IG-score is this value minus the DKL between the unigram frequency distribution in the predicted sequence and the corresponding distribution in high resolution protein structures from the PDB. This results in a single scalar DKL. The latter values are sourced from the publicly available trDesign code repository6.
5https://www.github.com/psipred/DMPfold2 6https://github.com/gjoni/trDesign
A.2 DARK0-Grad: SAMPLING SEQUENCES FROM DMPFOLD2 WITH GRADIENT DESCENT
We wished to investigate whether the seed sequences we use in DARK could be generated by a minimizing the negative IG-loss predicted by DMPfold2 with respect to the input sequence. We do so implementing and testing an approach similar to those shown by Killoran et al. (2017) and Norn et al. (2021), but instead using DMPfold2. We begin with an input of shape [B,L, 22] where the first dimension is the number of sequences in the minibatch. To keep compute time results comparable to the other methods presented here we use B = 100. The input itself is a one-hot encoded random sequence of L = 100. The version of DMPfold2 we use is altered to predict distograms for batches of sequences under the constraint that the alignment is always of size 1, which is not a constraint for us here.
To begin optimization, we compute the distograms from the batch of input sequences, backpropagate the negative IG-score as our loss, and then update the original input [B,L, 22] vector. This is one step. As we are back propagating through a discrete node, and the updated input vector is no longer discrete, we do the following in all optimization steps to make the operation differentiable. We softmax the updated input vector (we refer to this now as a profile) in the final dimension and then we sample sequences from the softmaxed residue distributions using a straight-through estimator (Bengio et al., 2013).
We perform 200 of the above steps as we found it sufficient for the negative IG-scores to converge. Instead of sampling final sequences from the final profiles, we track and keep the sampled sequence that achieves the highest IG-score during the 200 steps. We explored a variety of optimizers, profile regularization and normalization approaches, and other variants with the aim of improving the average final IG-score. Despite this, we generally found Adam (Kingma & Ba, 2014) with a learning rate of 1, default β parameters, and a temperature of 1.0 on the softmax to be the most consistently effective approach. To calculate the average IG-score of this method, and the average sample time, we generate 1000 sequences from the described approach. These resulting metrics can be seen in Table 1.
A.3 DARK1-Adversarial:
As part of our investigations we were motivated to assess the effect of incorporating a structure predictor directly into the training of our autoregressive DARK models. There has been a range of recent work (Larsen et al., 2016; Esser et al., 2021) exploring the combination of learned and fixed discriminators, like the learned descriminator in generative adversarial networks (GANs) (Goodfellow et al., 2014), with generative autoencoding models like variational autoencoders (VAEs) (Kingma & Welling, 2013; Rezende et al., 2014), vector quantized VAEs (VQ-VAEs) (Oord et al., 2017), and many others. Although this topic has mostly been explore with latent variable models, we wishes to investigate whether such an approach would benefit our autoregressive approach. We note however that there are potentially problems with this approach which we discuss below.
To explore this model, we began with the full training procedure and architecture of the DARK1 model. We then incorporated DMPfold2 by, during a forward pass in training, taking the output probabilities for the batch in question and sampling one-hot sequences combined with a straightthrough estimator like in Appendix A.2. The negative IG-score is then calculated as a loss, LIG. We see using DMPfold2 in this way as a kind of adversarial regularizer, thus we termed it DARK1Adversarial. It must be noted that the DMPfold2 parameters are not being updated in this procedure, the model is effectively frozen. We calculate the final loss Lfinal as,
Lfinal = LAR + γLIG (3) where LAR is the original autoregressive cross-entropy loss and γ is a tunable hyperparameter. We performed a small search of γ values and found γ = 0.5 to perform well however this was not rigorously optimized. As can be seen in Table 1, the approach does not perform significantly differently from the original DARK1 model. We in part suspect that this is due to the DMPfold2 model simply reinforcing the loss signal from the cross-entropy against the original sequence; the original sequence is already optimized to have a high IG-score. A less charitable view is that the gradient returning to the decoder model from DMPfold2 is simply contributing noise. As mentioned, there is also a potential issue with the method itself, being that by backpropagating through DMPfold2, which sees the full sampled sequence, we are potentially exposing the decoder to gradient information from the full sequence, breaking the autoregressive nature of the model. Given DMPfold2
is frozen, and we did not see model collapse during training, we do not believe this is occurring, however, we highlight it as necessitating further empirical investigation.

heading: A.4 AUTOREGRESSIVE MODELS
text: We build the autoregressive neural network models used in this work with PyTorch (Paszke et al., 2019). We also use several utilities from the NLTK (Bird et al., 2009) package. For all neural network architectures we do mini-batch gradient descent against negative log likelihood (NLL) loss. We use the Adam optimizer (Kingma & Ba, 2014) with default β values of 0.9 & 0.999. Early stopping is performed using the validation set. A batch size of 64 is used with a learning rate of 1 × 10−3 for the RNN baseline, and 1 × 10−5 for self-attention models. The learning rate was not changed according to a schedule as it was not found to provide significant benefit, and so was left static. All training was done with automatic mixed precision as it brought a significant improvement in computational efficiency. Below we describe the architectures for the best RNN baseline model and for the self-attention models.
Inputs For all language models trained on synthetic data, the input is a batch of B sequences in an array sized [B,L]. This contains 22 integers corresponding to 20 amino acids, a start token, and a stop token. The length L is the length of the protein sequence with a start token placed at the beginning of the sequence. As is standard for autoregressive models, the target sequence used for calculating the loss is the protein sequence with a stop token placed at the end. All neural networks use embedding layers that use these as integer indices for embedding vectors learned during training.
RNN architecture The best RNN model uses an embedding layer, learned during training, of size 64. The embedded sequence of size [B,L, 64] is fed into an LSTM model with a hidden size H of 128 and 1 layer (Hochreiter & Schmidhuber, 1997). The hidden vector output at each time step of the LSTM is then fed to a linear layer (H = 128)→ ReLU activation→ linear layer (H = 22). The resulting final array, sized [B,L, 22], is then log-softmaxed.
Xfmr architecture The Big and Small Xmfr models use a common architecture but with different numbers of parameters. We describe the architecture using the sizes of the Big model first. An embedding layer of size 384 is concatenated with a positional encoding of size 384 (Vaswani et al., 2017), resulting in an array of size [B,L, 784]. The positional encoding uses a maximum length of 101. In comparison, the Small model uses an embedding size of 64, so the input to the following decoder has 128 dimensions in the final axis. We use a standard decoder architecture for autoregressive multi-headed self-attention with causal masking (Vaswani et al., 2017; Radford et al., 2019). This includes 10% Dropout (Srivastava et al., 2014) and the use of Layer Normalization (Ba et al., 2016). We use the Fast Transformers library (Katharopoulos et al., 2020) for implementing models. The Big model uses 12 layers with H = 784 and 12 heads, and the Small model uses 4 layers with H = 128 and 4 heads. For both models the output of the decoder goes into a linear layer (H = 128) → ReLU activation→ linear layer (H = 22). The output is produced is also a [B,L, 22] array that is log-softmaxed.
For both RNN language models and self-attention language models, additional decoding strategies beyond simple autoregressive sampling like Top-K sampling (Welleck et al., 2019) were not found to significantly improve results. We mask out cysteine when sampling from DARK models, however we include it in the sizing of the input embeddings and output size to maintain consistency with other approaches. The output array of size 22 accounts for the 20 amino acids, a start, and a stop token. Further investigation of more complex decoding schemes, as well as the use of different positional encoding methods and types of attention, are potential avenues for future research.

heading: A.5 TRAINING AN UNCONDITIONAL MODEL WITH NATURAL SEQUENCES
text: Here we describe further architecture and training details of the ‘UF50’ model used to benchmark the DARK models produced in this work. This model provides a direct comparison of training with natural sequences and training with synthetic proteins from DARK. UF50 is an unconditional language model that uses the same training procedure and architecture as the DARK3 model, with some modifications, and trained on sequences from nature. Specifically, the model is trained on the representative sequences for each cluster in the UniRef50 database (The UniProt Consortium,
2021), yielding a training set of 50,106,395 sequences. UF50 provides a benchmark for an approach similar to the ProGEN model (Madani et al., 2020), if ProGEN was an unconditional model. ProGEN is conditioned on a variety of evolutionary and functional information for the task of protein engineering, making it not directly comparable to the work performed here.
We now describe the modifications used in the UF50 model. During training, missing amino acids and non-canonical amino acids are mapped to an additional ‘UNK’ token which does not contribute to the loss. Thus, the final vocabulary size is 24. Batches are dynamically padded to the length of the longest sequence, and padded positions are masked for the loss calculation. To reduce the GPU memory required for long sequences, any sequences longer than 100 amino acids (start/stop not included) are clipped to 100. For the target sequence, instead of placing a stop token at position 101, the amino acid sequence in position 101 from the original non-clipped sequence is used. In this way, due to the causal masking, this is exactly the same information available to the model as if it was not using clipped sequences, within the range of 101 amino acids. As all evaluated sequences are of length 100, the clipping has no effect on any conclusions drawn. To provide a level playing field, UF50 uses the same validation set as DARK models. We also use rejection sampling to sample 1000 sequences of length 100 for evaluation and comparison. Additional decoding strategies like Top-K sampling and placing a temperature factor on the next-token softmax were explored but not found to aid performance.

heading: A.6 FURTHER COMPUTATIONAL EFFICIENCY AND HARDWARE DETAILS
text: We report all timing estimations, unless stated otherwise, based on using an RTX 2080 Ti graphics card on the same single machine. We report all sample speed timings in seconds per example. Sample timing estimations were calculated with batches of 100 sequences, and calculated on the same single machine. The only example were this was not the case was for DARK0-Grad which was performed in batches of 10 due to memory constraints. Model training and sample refinement was performed with resources available at the time but these typically consisted of V100 graphics cards. As such, we report timing for training, refinement, and simulated annealing experiments as they would take on V100 graphics cards. The average time taken for refinement is 51.3 seconds.

heading: A.7 EXCLUSION OF CYSTEINE
text: Cysteines are a special case with regards to the other amino acids as they are primarily present in proteins to form disulfide bonds. Disulfide bonds are very strong bonds that form between two cysteine residues. They are, for the most part, present only in proteins that are secreted by the cell or exist on the outward facing surface of the cell membrane. This is due to the fact that the disulfide bond cannot form in the reducing environment of the cytosol of cells. We note that for all the statements, as always with biology, there are exceptions to the rules. In de novo design studies, the inclusion of cysteine is understandably usually framed as the inclusion of disulfide bonds. Marcos et al. (2018) and Bhardwaj et al. (2016) provide recent examples of de novo designing with the intent to use disulfide bonds or simply explore if their inclusion provides benefit. Disulfide bonds can be very useful for adding stability to a protein but complicate protein expression and depending on the aim of the design study, it may be disadvantageous or advantageous. From a structural perspective, a disulfide bond can very much ‘make or break’ a candidate so it requires consideration. Its for this reason that the inclusion of cysteine and disulfide bonds is typically addressed in a methodologically specific way, Bhardwaj et al. (2016) being a direct example of that. This is relevant in our case because it clearly indicates how the default stance is generally to not include cysteine and disulfide bonds unless there is a reason to. Ultimately, its exclusion does not make the problems we tackle in this work any easier and if it were to be included it would be unlikely to add any significant value.

heading: B SEQUENCE AND STRUCTURE ANALYSIS
text: B.1 VALIDATING THE VALIDATION AND TEST SET
Ideally each seed sample from the black-box oracle could be viewed as identical and independently distributed (i.i.d.). However, given the sometimes unexpected pathologies (e.g. adversarial examples) that can be present in neural networks, we directly measure the diversity of the training set
sequences as well as look for any direct overlap of the validation and test sets with the training set. As can be seen in Table 4, the training set sequences begin to form clusters at around 30% sequence identity.
As a method of comparing generative models between architectures and iterations of DARK, a test set and a validation set of 950 samples each are generated, using the simulated annealing approach of trDesign. We check that the two sets do not share extremely similar or identical sequence samples with the seed set, producing artificially high performance metrics. As such, validation and test set sequences are separately compared to all 15,000 sequences in the seed sample set. This is done using the MMseqs2 (Steinegger & Söding, 2017) search tool with a sensitivity parameter of 7 and one iteration. Given the clustering results of the seed samples (Table 4), we expect a high proportion of the sequences to share some small similarity to sequences in the seed sample set. This was found to be the case; from the validation and test set, 854 out of the 950 sequences in each set were found to have some similarity to sequences in the training set. The mean sequence identity between the similar sequences for both sets was low, being 31% for both. The highest sequence identity found in the validation set was 68%, and for the test set 60%. These results are as would be expected; they suggest that both sets are sampled from the same distribution as the seed set while not containing any near-identical copies. This is acceptable, given that this is an unsupervised task using synthetic samples, where the aim is to learn that distribution. We make the distinction as had natural sequences been used a more stringent test and validation splitting based on structure and evolutionary relationships would be required.
As part of investigating the effects of DARK iterations, we also performed this search process with the test set for the 500,000 examples of the third iteration’s training set. The number of examples in the test set that had any detected hits increased a small amount, from 854 (N = 1) to 922. The mean sequence identity of the hits increased slightly from 31% to 33% and the highest identity hit increased from 60% to 66%. Again, these values are acceptable, and even potentially a positive sign. One likely cause for the increase in the performance of models across iterations of DARK is that the area of sequence space covered by the training set is progressively growing in an unsupervised manner. Even from a pessimistic perspective, if the N = 3 model has better likelihood estimates on the test set because the previous iterations generated training samples closer in space to the test set, then it suggests that the iterations are indeed expanding and filling in unseen parts of the distribution.
B.2 THE STRICT STRUCTURE-BASED TEST SET AND DARK2-STRICT
We take two key approaches to providing a stringent measure of how well DARK generalizes. Both are enabled by first using AlphaFold to generate all-atom structure predictions, and pLDDT scores, for all 100,000 sequences contained in the DARK2 training set. This set covers the 15,000 initial seed sequences and the 85,000 sequences from the first iteration of DARK. We also assign each structure a topology in the form of a topology string. This describes the topology by the secondary structure elements going from the N-terminus to C-terminus. For example, if this had been done for the de novo designed ‘Top7’ protein (PDB ID: 1QYS) the topology string would be ‘ββαβαββ’ or ‘bbababb’. We also repeat this procedure for the 950 sequence validation set.
The first way we use this 100K set is to group sequences by their topologies and show that the first iteration of DARK generates sequences with structures unseen in the initial 15K seed set. The seed set contains 421 unique topologies and the 85K set contains 1011 unique topologies, with an intersection of 373 between the two. As such, there are 638 topologies in the 85K set that were unseen in the training set. This accounts for 1476 sequences, of which 225 sequences (15.2%) have a pLDDT≥ 70. This provides strong evidence that the DARK iterations are expanding the training sets to include sequences with unseen and, importantly, ordered structures.
The second way we use this set is to construct, and train a model on, a new training and test set split based on structure, as is considered the gold-standard in evaluating structure prediction methods (Chothia & Lesk, 1986; Söding & Remmert, 2011). We refer to the model trained and tested on this new split as DARK2-STRICT. We construct the test set by first randomly sampling 101 topologies, approximately 10%, from the 1059 unique topologies in the 100K set, ignoring the 100 most common topologies and uninteresting topologies, being a single helix (131 sequences) and only coil (3 sequences). This 101 topologies accounts for 413 sequences which are removed from the 100K set, and become the strict sequence test set (See Appendix B.2 for the test topologies and further
details). We also remove and discard 4 sequences from the 950 example validation set that have a strict test set topology. The original 950 example test set is not used here.
Next we use the MMseqs2 search tool (Steinegger & Söding, 2017) with the maximum sensitivity parameter (-s 8.5) to aggressively search for potential similarity between the strict test set sequences and sequences in the remaining 99587 training examples and 946 validation examples. We then discard any training and validation sequences that MMseqs2 returns as hit, meaning it is potentially similar to a test set sequence. This reduces the size of the training and validation sets to 87551 and 621 examples respectively. We then use the 87551 sequence training set, 621 sequence validation set, and 413 sequence strict test set to train and evaluate the DARK2-STRICT model.
We note that this approach to assigning structures to topologies is an aggressive one as it does not consider the different packing arrangements that the secondary structure elements can have. Thus, our measures of the number of topologies present in the set is likely a lower bound on the number of topologies that would be assigned by a more intricate approach like that done for CATH (Orengo et al., 1997). This may make the test set we have constructed more difficult however we wish to take a conservative approach and challenge the model, so we find this acceptable.
The 101 randomly selected topologies, and their sequence counts, making up the STRICT test set The ‘a’ is for alpha helix and ‘b’ is for beta sheet. The 101 topologies cover a total of 413 sequences: bbabbbabbbb (7), abaaaab (2), baaababab (1), aaabbbabb (1), babb (3), bbbbbbbbbbbab (6), bababaaa (1), abbaabbbb (2), bbbabbababb (3), bbabbbaabb (3), bbabbbababbab (1), bbbbbbababb (3), aaabbb (5), babbabbbba (2), abbbbbbbabbab (1), aabbaa (32), bbaababa (1), bbbbbbbaab (4), bbabbbababb (1), abaaabbaa (1), abbaaabbaa (1), bbabbabbbbbabbb (1), abbaabba (7), babbbabab (3), babbabbabb (19), bbbbabbaabba (1), aabaaabbaa (1), babbbbbabbbb (2), bbbaabbba (2), bababbbbbbb (2), bbabbabbabbb (6), babbbaab (2), bbbababbb (18), babbbabba (3), bbbaabaa (1), babaabb (14), abaabbabbab (1), bbababbbb (19), baaaab (9), bbaaaabb (1), abaabaabaabaa (1), bbaabab (4), abbbababbbbb (1), bbbbaabbbb (1), bbbbbbbbaab (1), baabbaab (2), bbababababab (1), bbabbabbabbab (2), aabbbab (4), bbaabbabbb (1), bbbbbbbabab (5), bbaababaab (1), bbabbbbab (29), abbbbabb (22), aaaabab (2), babbbbbabbb (4), bababbabb (6), babaabbbaab (1), bbabbabbbbbbb (1), bbbbabbbabbb (2), babbabbabab (1), bbabaabbb (2), baabbbbbb (5), aaabaabaab (1), bbbaabb (7), abbbbabbabb (1), bbaabbbbaab (1), ababaaa (1), babaabab (7), bbabbbbaa (3), aabbaabb (4), abbababb (5), abbbabbabbb (1), bbaabaabb (1), bbbbaab (7), baababb (8), aaabaabba (1), abababaa (1), bbbabbbbbbabbab (1), bbabbbbabbbb (4), abbaaaaa (3), abbbbababb (1), bbbbabbbaab (1), abaabbabba (1), bbbbabbbbabbb (2), bbabbbbabbbabb (1), babbbba (26), aaabbabbb (1), abbbabbbbbbb (1), babbbbabab (3), bbaabbabbba (1), bbbbbabbabbb (1), abbabbabbb (2), bbbaabbbabba (1), bbabbbbbbbba (1), babbababb (8), abb (5), bbaaabbbbba (1), bababaab (2), bbbbabbbbbbbb (2), babbbbabba (3)

heading: B.3 PREDICTING STRUCTURE WITH ALPHAFOLD
text: We use the publicly available version of AlphaFold7 for predicting structures and as part of a quality metric. In both cases, we use the standard AlphaFold inference pipeline, meaning we predict with the five ‘CASP14’ models, take the predicted structure with the highest pLDDT as the final result, and perform relaxation. We use adapted inference scripts to accommodate for predicting with only a single sequence and without using the template pipeline or the multiple sequence alignment (MSA) pipeline. This consisted of setting all models to use the same default configuration, as they all use variants on the default, and then setting all to not use templates. The MSA clusters are reduced to 1 sequence, being the single input sequence, and the deterministic flag is set to true. The number of recycles is set to 4. Following structure prediction, the standard relaxation is performed for the best final structure. We note that relaxation is not necessary to produce a final structure, thus it can be said that a trained language model, like those we present, with AlphaFold provides a direct way to sample from a joint distribution of sequence and structure, and without additional conditioning information. Predicting structure with AlphaFold provides the pLDDT scores at each
7https://github.com/deepmind/alphafold
residue position. We take the average of these for the pLDDT score of a given sequence. Structures are visualized with ChimeraX (Pettersen et al., 2021).
B.4 ALPHAFOLD PREDICTION CONFIDENCE AND DARK3 FOR RECENTLY RELEASED DE NOVO DESIGNED STRUCTURES IN THE PDB
As we use the pLDDT from AlphaFold as a metric for measuring structural order, we wish to doublecheck that AlphaFold assigns good or better pLDDT scores to de novo designed sequences in the PDB that were not used as a part of its training set. We manually curate a list of 26 of these sequences from the PDB. The PDB IDs containing these are 6MRR, 6MRS, 6NUK, 6MSP, 6XWI, 6XXV, 6VTW, 7KBQ, 6D0T, 7JH6, 6XEH, 6X9Z, 6WVS, 6E5C, 6DS9, 5TX8, 5W9F, 6YWD, 6Z35, 6W6X, 6W90, 6WI5, 6W3W, 6VGB, 6W70, 6YQY.
This is not a comprehensive list as, amongst other reasons, some sequences were not included if they were very sequence-similar to an already selected structure or if they are a large tandem repeat structure (e.g. ‘armadillo repeats’) to avoid biasing the average. This is still a biased set of structures as they are the product of whichever goal a given study was investigating, but the pLDDT scores are still informative. We note that although the PDB has a ‘DE NOVO PROTEIN’ classification, this includes engineered, designed, and de novo designed proteins. The average calculated pLDDT of all sequences in the list is 88.1 (σ = 4.41), a good score, as is expected.
We also calculate the likelihood of these sequences under DARK3 compared to random sequences of length 100. We ignore the 5 sequences that contain cysteine. To account for DARK having been trained on static length proteins we ignore the prediction of the ‘stop token’ and ignore all residues in a sequence above position 100 if they exist. We calculate the bits per character (BPC), where characters are residues, across 100 random sequences and the 21 designed sequences. The random sequences have a mean BPC of 5.88 bits and the designed sequences have a mean BPC of 4.96 bits. In having a lower BPC suggests that the DARK models have learned a distribution of sequences that incorporates structure and is capable of producing de novo designed sequences more so than random. It is also unlikely that by random chance the DARK models have generated and learned from sequences in the same parts of the dark sequence space

heading: B.5 BENCHMARKING WEAKLY CONDITIONED METHODS
text: We provide results for benchmarking with Good+ pLDDT proportions on two recent weakly conditioned generative protein design models, which we refer to as gramVAE (Greener et al., 2018) and gcWGAN (Karimi et al., 2020). These both use a fold level description that we call weak as they are far removed from the explicit structural information provided in fixed backbone models (Ingraham et al., 2019; Strokach et al., 2020; Cao et al., 2021). We stress that these models are not directly comparable to the unconditional models explored in this work. Furthermore, we provide these only to gain a small degree of insight but this is significantly muddied by both models being able to ‘cheat’. If either of them reproduce a sequence from their training set, then it is extremely likely that the same sequence and its structure are in the AlphaFold training set8 and so likely to be assigned a higher pLDDT than would otherwise be expected. Nonetheless, we provide these benchmarks as they provide some small degree of context. Even with extra information from conditioning variables and this ability to cheat, they are still outperformed by DARK3 (Table 8). We speculative that this difference can even be seen visually in the predicted structures shown by both models (Figure 5c of Karimi et al. (2020) and Figure 6 of Greener et al. (2018)) which contain clear disordered regions.
The gramVAE9 uses a context-free grammar of SCOP protein folds (Andreeva et al., 2020) as a conditioning variable for a Variational AutoEncoder (VAE). Both encoder and decoder use nonlinear layers and sequences are padded to a maximum length of 140. The gcWGAN is a conditional Wasserstein Generative Adversarial Network that also pads length to use fixed-sized networks (Karimi et al., 2020). It conditions on a small latent representation of folds, learned from performing kernel PCA on an all-vs-all distance matrix between all folds in the dataset, where the kernel used is symmetric TM-score (Zhang & Skolnick, 2005) between two folds. The publicly available
8AlphaFold makes the blanket statement that all structures from the PDB with a maximum release date of 30 April 2018 are used in its training.
9https://github.com/psipred/protein-vae
implementation10 provides 600 sequences generated for the 6 folds (100 sequence each) in its test set using the gcWGAN and the gramVAE in a head to head comparison. Being publicly available and a direct comparison, we use these provided sequences to calculate Good+ pLDDT proportions for both models. The gramVAE achieves 21.0% and the gcWGAN achieves 2.83%.
We note that the majority of the folds in the gcWGAN test set are in the gramVAE training set, which likely explains part of the performance difference. We are not concerned by this as the DARK models still outperform the gramVAE with an artificial boost. We also note that the gcWGAN model does not discretize the output of the generator and use any approaches to backpropagating through the discrete output variable (Bengio et al., 2013; Jang et al., 2016). This provides the discriminator a trivial learning scheme: place the decision boundaries at the edge of the probability simplex of tokens at any position. This may explain part of the gcWGAN’s poor performance. In both cases, we recognize that these models are not ideally suited for sequence problems and have reported poor performance elsewhere (Cao et al., 2021). We expect future weak conditioning methods will provide more competitive benchmarks. Ideally, we hope that future investigations into novel unconditional models will provide direct benchmarks, establishing the best approach to unconditional modelling.

heading: B.6 DARK MODEL SAMPLES ARE DISTINCT FROM NATURAL PROTEINS
text: We find that samples from the DARK3 model are distinct from sequences in nature. A sequence search against UniRef100 of the 1000 evaluated sequence samples (without refinement) returned no hits for all sequences except one, which had 5 weak hits. Two pairs of these 5 sequences were extremely sequence-similar. Repeating the sequence search with a different tool, HHblits (Remmert et al., 2012), no natural sequences were found that met the HHblits default threshold value (an E-value below 1× 10−3).

heading: B.7 SEQUENCE SEARCHING AND CLUSTERING
text: We use the MMseqs2 package (Steinegger & Söding, 2017) to perform sequence searches and sequence clustering. For clustering, the default clustering pipeline and settings are used with sequence identity as the threshold. If not specified, this is assumed to be a threshold of 30% sequence identity. For sequence searching we use the default pipeline of MMseqs2 and a high sensitivity parameter of 7.0 with 1 iteration against UniRef100 (03-2021 release) (The UniProt Consortium, 2021). Any sequence searching done with HHblits (Remmert et al., 2012) uses an E-value cut-off of 1 × 10−3 with 1 iteration against UniRef30 (06-2021 release) (Mirdita et al., 2017). In this work, we call any two sequences similar if MMseqs2 (or HHblits) returns a hit between them, meaning it has detected similarity. Sequences that have a, for example, 40% sequence identity or above are also defacto assumed as similar.

heading: B.8 STRUCTURE CLUSTERING
text: To cluster predicted protein structures we perform single-linkage agglomerative clustering using TM-score (Zhang & Skolnick, 2005) as a similarity metric, meaning an all-against-all matrix of TM-scores is calculated for each structure in the set of 1000 samples. We use a TM-score of 0.4 as the cutoff for selecting best cluster centroids and a TM-score of 0.65 for single-linkage clustering around the centroids.

heading: B.9 GENERATING DE NOVO DESIGNS WITH DARK AND ALPHAFOLD REFINEMENT
text: Background To give a concrete demonstration of how DARK3, and models like it, can be used in developing novel methods, we develop a simple but novel design method. The method we develop, termed AlphaFold refinement, tackles what could be considered the next step for using a model like DARK. Although we have used AlphaFold primarily to evaluate samples, and by extension models, it can also be used to improve candidates. An obvious first step is that AlphaFold produces both distogram based predictions and an all-atom structure prediction, so a simulated annealing approach like trDesign’s (Anishchenko et al., 2020) could be used to sample sequences. We note that the
10https://github.com/Shen-Lab/gcWGAN
all-atom prediction can be considered a point estimate of the distogram, but the all-atom model is ultimately the output of the most interest to us for design.
The challenges of a simulated annealing approach with AlphaFold Unfortunately, this approach poses two problems. The first is that, at its fastest, meaning running AlphaFold with no recycling (Jumper et al., 2021), each forward pass takes approximately 5 seconds. This is in contrast to 0.064 seconds with DMPfold2. Performing a simulated annealing run with the same schedule as we use would take ∼2.3 days which is prohibitively long. The second problem is one we find in extending the approach itself. Directly optimizing on the AlphaFold IG-score (AIG-score) itself appears to produce, what amounts to adversarial attacks on AlphaFold. This can be seen in Figure 2A & 4, and is particularly evident in the low AIG-scores (low compared to our approach) but extremely high pLDDT scores. The visual confirmation of this is evident in Figure 4 where all structures produced are a very unrealistic single long helix. We implement the AIG-score as the IG-score but using AlphaFold to calculate the distogram prediction from which the p(y|x) term is calculated. As we do not have the resources that were use to train AlphaFold (Jumper et al., 2021), we calculate a proxy background distribution by predicting the structure of 1000 random sequences of L = 100 and then averaging their resulting distogram predictions together into one vector.
Using DARK3 to enable AlphaFold design Introducing DARK3 into this problem brings several benefits. First and foremost it allows us to generate diverse and high confidence samples from AlphaFold without them being adversarial attacks. It also reduces the number of design steps such an approach would take, and adds a measure of control into the design process. We call this method AlphaFold refinement as it is similar in concept to the refinement done with DMPfold2 between iterations of DARK. DARK3 allows us to generate large sequence databases, with predicted structures, that have a diverse range of topologies. These include structures with topologies that have yet to be successfully de novo designed, and many with high confidence stable structures that can be used as de novo scaffolds into which functional sites are grafted (Woolfson, 2021). In short, DARK3 produces sequences that have structures that are of design interest.
AlphaFold refinment AlphaFold refinement begins with a sequence sample from DARK and its AlphaFold predicted structure, chosen by the designer. This is how control is introduced into the design process. For a graphic displaying the full AlphaFold refinement procedure see Figure 2. The next step is to take the predicted all-atom structure and calculate its one-hot encoded distogram representation (Jumper et al., 2021). We then bootstrap this representation as a fixed optimization target. The actual optimization that occurs is a greedy hill-climb on the sequence against the CrossEntropy between the predicted distogram and the bootstrapped target. As we use all five parameter sets when we use AlphaFold, we treat the distogram prediction at each step as if it were an ensemble of five AlphaFolds; we average the five predicted distograms as one final predicted distogram with which to calculate the Cross-entropy. We run this for AlphaFold refinement for 10000 steps, as can be seen in Figure 2C, however we found comparable results could be achieved in 7500 steps, being 10 hours per sequence on a single consumer-grade GPU (RTX 2080 Ti). This is significantly more efficient than an approach like using Rosetta which can take thousands of CPU hours. In summary, we mutate the sequence and keep mutations that minimize the Cross-Entropy to the original structure. The underlying reason for undertaking this approach is that the structure and sequence we are interested in will, more often than not, have a lower pLDDT than is possible with a similar sequence that folds to the same structure. In short, there is always room for improvement, and AlphaFold with DARK3 provides a way to address this.
This provides two benefits. The first is that it forces the objective to keep the original structure we are interested in. We find that this also provides a shorter and, we speculate, more direct optimization path to a high confidence structure. There stands to reason as we giving a fixed objective, a much easier task, rather than letting the random walk find arbitrary minima. This is evident in the average AIG-score and pLDDT score improving vastly while the final sequences remain very similar from a sequence identity perspective. What is particularly notable is that the average pLDDT score of the 10 examples explored is very similar to that of recently released de novo designed structures in the PDB (Figure 2D & 2E). The 10 examples we take are chosen from the 5000 DARK3 samples previously generated in a semi-random we fashion. We chose them randomly but rejected any with the exact same topology, accounting for packing, and any single sheet (CATH 2.20; Orengo et al. (1997)) topologies simply because they are uninteresting. The 10 examples cover 4 α& β structures,
3 α structures, and 3 β structures. This rather even spread of examples was entirely coincidental but convenient.
Another aspect of using the distogram is that it provides a softer means of keeping the original structure. In essence, by working in a probabilistic space, AlphaFold can perform small rearrangements to the original structure to improve its confidence while not straying to far. This is difficult to quantify but can be seen in Figure 2D Example 8, a complex α & β topology where several of the sheets are brought together in the lower half of the structure (as visualized) to improve packing. To summarize, in the limit of AlphaFold’s prodigious accuracy, AlphaFold refinement provides a way to refine design candidates of interest and improve the confidence (and order). Another way of viewing this approach is as an entirely neural network driven approach flexible-backbone design, here used as part of an overall de novo design method. In that vein, we note that this approach could be extended very easily to other tasks in protein design like fixed-backbone design.
Designing a novel immunoglobin As an example of using this approach we take a immunoglobin from the 5000 DARK3 samples that we previously generated and then AlphaFold refine it. We take the first example of an immunoglobin in the set and find that we successfully refine its structure to a high pLDDT and AIG-score while not altering the original approach. The results of this design can be seen in Figure 3. This provides a concrete example of using DARK3 and AlphaFold to design a high confidence candidate for a topology of high scientific interest that traditional approaches like Rosetta are known to struggle with.

heading: DARK Model Good+ pLDDT (%) ↑ Sequence Clusters ↑ Structure Clusters ↑

heading: C ADDITIONAL FIGURES AND TABLES

heading: Language Model Good+ pLDDT (%) ↑ Conditioning Training Data

Review summaries:
I think it helps to start by summarizing [Anischenko et al. 2020, Protein design by deep network hallucination], and then explain the differences in the submission.  [Anischenko et al. 2020] aim to generate AA sequences that are diverse but that still reflect the properties of training structures. * Their goal is to generate protein sequences whose structures tend to be unlike those of random sequences. Since random sequences have structures that tend to be unordered, they incidentally generate sequences having ordered structures. * They use the deep resnet from rtRosetta (similar to AlphaFold's) to predict interresidue geometry maps from sequence, and that is how structure is represented and compared. * Their design goal is to design samples from p(structure | sequence) p(sequence). * To represent p(structure | sequence) they use a proxy objective: the KL divergence between the predicted structures and the structures of random sequences; in other words, good sequences are those whose distance maps are very different from the distance maps of random sequences. * To represent p(sequence) they just use -KL(AA frequences of predicted sequences || AA frequences of PDB sequences). * Sequences are designed by starting from random and then annealing according to the proxy objective. * Designed sequences are then clustered by template modeling score (TM-score) and selected (from each cluster) based on  a measure of prediction 'consistency'.   Now, the submitted manuscript (DARK) is the same objective and proxy objective as [Anischenko et al. 2020], but with following methodological changes: 1. For the sequence->structure oracle, instead of using rtRosetta's interresidue geometry maps as the structure representation, here distograms [Senior et al. 2020] are used, and are output from a modified version of DMPfold2. 2. Rather than relying solely on one round of "random initialization -> annealing", DARK uses an iterative self-training scheme. First they do one pass of [Anischenko et al. 2020] to get 1st round of annealed sequences. The annealed sequences from 1st round are used to train an autoregressive generative model (Transformer-based). That model is sampled from to generate additional "random initializations" (which are now guided, not purely random). These guided samples are then optimized (just hill-climbing now instead of annealing, for speed), and added to the training set for subsequent rounds, etc.  There are other differences from [Anischenko et al. 2020], like using pLDDT score (orderedness) as a measure of "success", but the above two seem to be the main differences. The first of these two is worth trying but is not a conceptual contribution. The main idea of the paper is the second: self-training to generate diversified 'intermediate' initializations for further refinement. Basically it's using the initial refined sequences as a self-prior for exploring new sequences to refine.

The paper present a method to sample jointly de novo protein sequences and structures without relying on any structure conditioning. The main contribution is a new iterative procedure to generate a dataset of increasing size in an unsupervised fashion. This dataset is finally used to learn a generative model that captures the distribution over the "unnatural" sequence space, structures are predicted by using AlphaFold on the generated sequences. The authors provide empirical evidence of the quality (diversity and structural stability) of the de novo designs.

The authors propose a way to generate novel protein sequences that are likely to have folded structures. To do so, they start with a pre-trained sequence-to-structure predictor, and use simulated annealing to generate a set of initial sequences that are likely to have a structure. These seed sequences are used to train an unsupervised sequence model (an autoregressive transformer). More sequences are generated from this sequence model, and they are optimized to maximize the likelihood they have structures (using the pre-trained sequence-to-structure model), and are added to the growing set of training sequences. This process is repeated for several iterations to finetune the unsupervised sequence model. The authors show that sequences generated from this sequence model are likely to be structured as predicted by AlphaFold.