Year: 2022

Abstract:
Black-box optimization formulations for biological sequence design have drawn recent attention due to their promising potential impact on the pharmaceutical industry. In this work, we propose to unify two seemingly distinct worlds: likelihood-free inference and black-box optimization, under one probabilistic framework. In tandem, we provide a recipe for constructing various sequence design methods based on this framework. We show how previous optimization approaches can be “reinvented” in our framework, and further propose new probabilistic black-box optimization algorithms. Extensive experiments on sequence design application illustrate the benefits of the proposed methodology.

heading: 1 INTRODUCTION
text: Discovering new drugs to fulfill specific criteria, such as binding affinity towards a given molecular target, is a fundamental problem in chemistry and the pharmaceutical industry (Hughes et al., 2011). In this work, we focus on an important subdomain: de novo biological sequence design. This task is challenging for two reasons: (1) the exploration space for sequences is combinatorially large; and (2) sequence usefulness is evaluated via a complicated process which usually involves time-consuming and expensive wet-lab experiments.
Despite the difficulty of this task, many approaches have been developed over the past few decades thanks to recent advances in biochemistry and machine learning. The Nobel Prize wining paradigm, directed evolution (Chen & Arnold, 1991), which conducts local evolutionary search under human guidance, is one of the popular techniques. Unfortunately, it is limited by its sample inefficiency and reliance on strong prior knowledge, e.g., about where to mutate (Ahn et al., 2020). Furthermore, to compete with other machine learning methods (Gottipati et al., 2020), guided evolution (Yoshikawa et al., 2018; Jensen, 2019; Nigam et al., 2019) heavily relies on human intuition for designing domain-specific evolutionary operators, which may not always apply to tasks at hand.
In this work, we deem sequence design to be a black-box optimization problem, tasked with maximizing an unknown oracle function. We assume that oracle queries are limited due to the constraint on resources, such as the budgets for evaluating queries in a wet-lab. Thus, sample efficiency is crucial. We develop a probabilistic framework by reformulating the aforementioned black-box optimization target as a posterior modeling problem. With this framework, we draw a surprising connection between likelihood-free inference and sequence design, and thus linking two fields which are previously considered as unrelated. The key observation we leverage here for establishing this connection is that both settings share similar elements and targets which will be elaborated in Section 2.2. This connection facilitates our understanding of both fields and provides a recipe for developing sequence design algorithms. Going beyond, we also combine different probabilistic modeling insights and develop three novel composite probabilistic algorithms. We point out that our framework could actually be applied to any black-box optimization settings, but in this work we focus on its application to biological sequence design.
To demonstrate the empirical effectiveness of our methods, we conduct systematical experiments to evaluate their performance on four in-silico sequence design benchmarks. Our proposed meth-
∗Corresponding Author
ods achieve at least comparable results to existing baselines, and the proposed composite methods behave consistently better than all other ones across various sequence design tasks.
We summarize our contribution as follows:
• We develop a probabilistic framework that unifies likelihood-free inference and black-box optimization.
• Based on this framework, we provide a recipe for designing algorithms for black-box problems. We apply these ideas to propose a series of composite design algorithms.
• We perform systematical evaluation on a series of black-box sequence design benchmarks, and find that these algorithms achieves consistently comparable or better results compared to previous ones, thus illustrating the benefit of the proposed unified framework.

heading: 2 A UNIFYING PROBABILISTIC FRAMEWORK

heading: 2.1 BACKGROUND
text: Likelihood-free inference (LFI). We use θ ∈ Θ and x ∈ X to separately denote the parameters and the data generated via the mechanism x ∼ p(x|θ). In this scenario, LFI refers to a special kind of Bayesian inference setting where the likelihood function is not tractable but sampling (by simulation) from the likelihood is feasible. Consider the objective of modeling the Bayesian posterior when we cannot compute the likelihood p(xo|θ):
p(θ|xo) ∝ p(θ) p(xo|θ)︸ ︷︷ ︸ ? , (1)
where xo is the observed data, p(θ) is the (given) prior over the model parameters θ, p(x|θ) is the intractable likelihood function and p(θ|x) is the desired posterior over θ. While we do not have access to the exact likelihood, we can still simulate (sample) data x from the model simulator: x ∼ p(x|θ). Instead of trying to obtain a numerical value of the generic posterior p(θ|x) for arbitrary x, LFI only tries to obtain an approximation of p(θ|xo) for the given xo. During the inference process, we can take advantage of the sampled data: D = {(θi,xi)}ni=1 where xi ∼ p(x|θi) for selected values of θi.
Biological black-box sequence design. We consider biological sequence design as a black-box optimization problem:
m∗ = arg max m∈M f(m),
where f(·) is the oracle score function, and we would like to discover values of m for which f(m) is large. In real-world situations, a query of this oracle f could represent a series of wet-lab experiments to measure specific chemical properties or specificity for a given binding site target. In general, these experiments are time- and cost-consuming. As a result, the total number of queries is limited.
In our setting, we useM = VL to denote the search space for sequences with fixed length L, where V is the vocabulary for each entry of the sequence: for DNA nucleotides |V| = 4, and for protein amino acids |V| = 20. For variable length setting, we haveM = ∪L∈[Lmin,Lmax]VL, where Lmin and Lmax are the minimal and maximal length, respectively.

heading: 2.2 CONNECTING LFI AND BLACK-BOX OPTIMIZATION
text: In order to draw a connection to LFI, we require a probabilistic formulation of the black-box sequence design problem. To this end, we relax the goal of searching for a single maximum of the oracle / score function f to a posterior modeling problem, i.e., finding a representative sample of the configurations of m sampled with probability related to some target posterior. Think of C is the set of sequences with these desirable configurations, E is a Boolean event about whether a sequence m belongs to C, and our goal is to characterize the posterior distribution p(m|E) from which we obtain the desired sequences. Below, we consider two specific ways of doing this:
Example A. We explicitly define C (and E accordingly) as all the sequences whose scores are larger than a given threshold s:
C = {m|f(m) ≥ s}. (2)
Here s could be any fixed value, or a certain quantile of a particular score distribution. In this way, we have p(E|m) = p(m ∈ C|m) = 1{f(m) ≥ s} where 1{} is the indicator function. Example B. In a softer version of E and C, we can define its conditional probability of being true to follow a Boltzmann distribution:
p(E|m) = p(m ∈ C|m) ∝ exp(f(m)/τ). (3)
where τ is a temperature parameter. We introduce the exponential because f(·) does not necessarily take positive values. Any monotone transformation of f(·) to non-negative reals could be used, so that sequences with larger oracle scores have a greater probability of making E true.
With this posterior objective, our goal now becomes effectively modeling and sampling from the posterior p(m|E). It is thus natural to resort to the tools of Bayesian inference for this task. In order to examine this possibility, we draw a detailed comparison between the settings of black-box sequence design problem and likelihood-free Bayesian inference in Table 1.
It can be observed that both tasks share similar elements and targets. The two settings also share similar limitations on the allowed queries, which are too time-consuming and / or cost-intensive. Notice that in sequence design, the oracle could be either exact or noisy, thus we use the more general s ∼ f(m) formulation rather than s = f(m). We will further present several concrete examples as demonstrations of this correspondence in the following section.
Another way to understand this correspondence is to consider the following mapping T :
T : Θ×X →M× R (θ,x) 7→ (m, s), s.t. s = −‖x− xo‖ .
Here we can see the score value s as a quantitative metric for how close the generated data x (given θ) is to the target observed data xo. In addition, querying the oracle in the sequence design setting can also be thought of as follows: (1) sample x ∼ p(·|θ) and then (2) calculate s = −‖x− xo‖ under some distance ‖·‖. In this manner, T could conceptually transform any LFI problem into a black-box optimization task. In this work, we only focus on the application of sequence design.

heading: 3 METHODOLOGY
text: We provide a recipe for designing new sequence design algorithms based on the correspondence in Section 2.2. The recipe induces different approaches by modeling different probabilistic component of the Bayesian inference problem. We begin with common algorithm restrictions under this setting.
Common constraint for algorithms. Due to the restriction of simulation / query in our setting, we constrain our algorithms to act in a sequential / iterative way, gradually achieving the desired posterior round by round. Every algorithm starts with an empty dataset D = ∅ and an initial proposal p1(·) = p(·), where p(·) is the prior given by the task. In the r-th round of this multi-round setting, the algorithm would use the proposal pr(·) of this round to sample a batch of data (θ / m) for simulation / query, and augment the current datasetD with the newly obtained batch of data. We use n to denote the batch size for each round’s simulation / query. Afterwards, the algorithm updates
the proposal to pr+1(·). The outcomes for the two settings we discuss may be slightly different: an algorithm for likelihood-free inference would return the posterior, while a sequence design method would return the dataset of all the sequences it has queried, which hopefully contains desired high scored sequences. On the other hand, a sequence design method could produce as an intermediate result a generative model for sampling queries, which then completely fits with the LFI framework.

heading: 3.1 BACKWARD MODELING OF THE MECHANISM
text: Approximate Bayesian Computation (ABC) (Beaumont et al., 2002) is a standard method for tackling LFI problems. In Algorithm 1, we display one of the most popular variants: Sequential Monte Carlo-Approximate Bayesian Computation (SMC-ABC) (Beaumont et al., 2009). In each round, parameters θ are sampled from the current proposal distribution pr(θ) for simulation. A rejection step is then involved to remove the θi whose simulation outcomes xi cannot reproduce the observed data xo with sufficient accuracy. The remaining accepted {θi}i are adopted to update the next round’s proposal pr+1(·) towards the target posterior, i.e., by refitting qφ with the modified data. We defer more details of this approach to Section A.1 in Appendix.
Algorithm 1 SMC-ABC p1(θ)← p(θ); for r in 1 to R do
repeat sample θi ∼ pr(θ); simulate xi ∼ p(x|θi); until n samples are obtained D ← D ∪ {(θi,xi)}ni=1 sort D according to −‖xi − xo‖; fit qφ(θ) with top {θi}i in D; pr+1(θ)← qφ(θ);
end for return p̂(θ|xo) = pR+1(θ)
Algorithm 2 FB-VAE p1(m)← p(m); for r in 1 to R do
repeat sample mi ∼ pr(m); query the oracle: si ← f(mi); until n samples are obtained D ← D ∪ {(mi, si)}ni=1; sort D according to si fit qφ(m) with top {mi}i in D; pr+1(m)← qφ(m);
end for return {m : (m, s) ∈ D}
It would then be natural to construct an analogical sequence design algorithm using top scored entities {mi}i to guide the update of a certain sequence distribution, see Algorithm 2. Interestingly, this is the proposed sequence design algorithm in Gupta & Zou (2019), where the authors name this kind of updating “feedback” because training of the parametric generator qφ(m) exploits feedback signals from the oracle. In this paper, we follow Brookes & Listgarten (2018) to crystallize qφ(m) to be a variational autoencoder (Kingma & Welling, 2014), and use the term Feedback-Variational AutoEncoder (FB-VAE) to refer to Algorithm 2. We place Algorithm 1 & 2 side-by-side to highlight their correspondence. We also make the same arrangement for the following Algorithm 3 & 4, Algorithm 5 & 6 and Algorithm 7 & 8.
Algorithm 3 Sequential Neural Posterior p1(θ)← p(θ); for r in 1 to R do
repeat sample θi ∼ pr(θ); simulate xi ∼ p(x|θi); until n samples are obtained D ← D ∪ {(θi,xi)}ni=1; qφ ← arg minq Ex [DKL(p(θ|x)‖q)]; pr+1(θ)← qφ(θ|xo);
end for return p̂(θ|xo) = pR+1(θ)
Algorithm 4 Design by Adaptive Sampling p1(m)← p(m); for r in 1 to R do
repeat sample mi ∼ pr(m); query the oracle: si ← f(mi); until n samples are obtained D ← D ∪ {(mi, si)}ni=1; qφ ← arg minqDKL(p(m|E)‖q); pr+1(m)← qφ(m);
end for return {m : (m, s) ∈ D}
In comparison with SMC-ABC, the Sequential Neural Posterior (SNP) method (Papamakarios & Murray, 2016; Lueckmann et al., 2017; Greenberg et al., 2019) for likelihood-free inference adopts a more flexible approach, taking the power of conditional neural density estimator (e.g., Papamakarios et al. (2017)) to model the general posterior p(θ|x), which takes arbitrary θ and x as two inputs and outputs a distribution. This neural estimator is trained via approximately minimizing the
Kullback–Leibler (KL) divergence between qφ(θ|x) and the true posterior p(θ|x). We defer more training details to Section A.1 in Appendix. Under the connection viewpoint, one similar algorithm for sequence design is the Design by Adaptive Sampling (DbAS) proposed in Brookes & Listgarten (2018) which is characterized in Algorithm 4, fitting qφ(m) through minimizing the KL divergence with the posterior p(m|E). Based on the difference in specific implementations, both algorithms have more than one variant, whose details are deferred to Section A.1 in Appendix.
We refer to the above algorithms as “backward modeling” because the trained generative network qφ (going from x/E to θ/m) is a sort of reverse model of the simulation mechanism (which goes from θ/m to x/s).

heading: 3.2 FORWARD MODELING OF THE MECHANISM
text: Whereas the above methods focus on directly modeling the target posterior with a generative model that learns a “reverse mechanism” of the simulation process, it is also possible to model the “forward mechanism”, which is consistent with the simulation process. Papamakarios et al. (2019) claim that the forward modeling approach may be an easier task than its backward counterpart, as unbiased estimation of the likelihood does not depend on the choice of proposal. Consequently, in contrast to SNP, Papamakarios et al. (2019) chooses to train a neural density estimator to model the conditional likelihood distribution qφ(x|θ) sequentially in each round. The training is achieved by maximizing the total log likelihood maxq ∑ i log qφ(xi|θi) with data samples from the dataset D at the current (r-th) round. The downside of this forward approach is an additional computational Markov Chain Monte Carlo (MCMC) step is needed to sample from the r-th round posterior / proposal pr(θ). The resulting approach, which is coined (Papamakarios et al., 2019) the Sequential Neural Likelihood (SNL), is summarized in Algorithm 5.
In the spirit of directly modeling the forward mechanism of sequence design, we train a regressor f̂φ(m) in a supervised manner to fit the oracle scorer. In order to adapt this regressor into the update procedure of the proposal of the next round, we use q̃(m) to denote the unknown posterior p(m|E) with knowledge of f̂φ(m) and prior p(m). The specific construction of q̃(m) depends on the choice of E . For instance, if we choose Example B in Section 2.2 to be the definition of E , then q̃(m) is the distribution with (unnormalized) probability p(m)·exp(f̂φ(m)/τ). See Section A.2 in Appendix for more elaboration about this point. We then choose the update procedure of the proposal pr+1(m) to be analogical to that of SNL. We name this proposed algorithm to be Iterative Scoring (IS) to avoid confusion with likelihood-free inference algorithms. Furthermore, depending on different definition of E , we use the name “IS-A” and “IS-B” for them in the following sections.
Algorithm 5 Sequential Neural Likelihood p1(θ)← p(θ); for r in 1 to R do
repeat sample θi ∼ pr(θ); simulate xi ∼ p(x|θi); until n samples are obtained D ← D ∪ {(θi,xi)}ni=1 fit qφ(x|θ) with D; pr+1(θ) ∝ p(θ) · qφ(xo|θ);
end for return p̂(θ|xo) = pR+1(θ)
Algorithm 6 Iterative Scoring p1(m)← p(m); for r in 1 to R do
repeat sample mi ∼ pr(m); query the oracle: si ← f(mi); until n samples are obtained D ← D ∪ {(mi, si)}ni=1; fit f̂φ(m) with D; construct q̃(m) with f̂φ(·) and p(m); pr+1(m)← q̃(m);
end for return {m : (m, s) ∈ D}

heading: 3.3 MODELING A PROBABILITY RATIO
text: In this subsection, we discuss yet another approach, through the estimation of a probability ratio. Gutmann & Hyvärinen (2010) proposes noise contrastive estimation as a statistical inference approach. This methodology turns a hard probability modeling problem into binary classification, which is considered easier to learn. In contrast to the aforementioned likelihood-free inference methods which rely on a form of density estimation to perform the task, Sequential Neural Ratio
(SNR) (Hermans et al., 2019) takes a similar approach as noise contrastive estimation. SNR adopts a classification approach to estimate the likelihood-to-evidence ratio r(θ,x) = p(θ|x)/p(θ). SNR is summarized in Algorithm 7. Specifically, in each round, SNR fits a binary classifier dφ(θ,x) ∈ [0, 1] in the following manner:
arg min d {∑ (θi,xi)∈D [− log d(xi,θi)] + ∑ (θ′i,xi)∈D′ [− log(1− d(xi,θ′i))] } . (4)
We show that with the D and D′ established in Algorithm 7, we have
d∗(θ,x) = p(θ|x)
p(θ) + p(θ|x) , r∗(θ,x) :=
d∗(θ,x) 1− d∗(θ,x) = p(θ|x) p(θ)
where the “∗” denotes the optimality. We have the following Proposition 1: Proposition 1. Let p0(a) and p1(a) be two distributions for da-dimension random variable a which takes value in the space of A = Rda , and d(a) : A → [0, 1] is a real-value function mapping any a to a positive real value number. Then the functional optimization problem
arg max d:A→[0,1]
{ Ea∼p0(a)[log d(a)] + Ea∼p1(a)[log(1− d(a))] } .
will lead to the optimal solution d∗(a) = p0(a)p0(a)+p1(a) .
See the proof in Section A.3 in Appendix. After training d, SNR can obtain the posterior density value by p(θ|x) = r∗(θ,x)p(θ). SNR mirrors SNL in that it samples the new proposal pr+1(θ) without explicitly modeling the posterior.
On the other hand, we propose a sequence design algorithm analogous to SNR and named Iterative Ratio (IR), which estimates the probability ratio r(m) = p(m|E)/p(m) for the purpose of posterior sampling. IR first builds two datasets D and D′, corresponding to two different distributions p(m|E) and p(m). We take a similar binary classification approach, whose training objective is mind {∑ m∈D[− log d(m)] + ∑ m∈D′ [− log(1− d(m))] } . The desired ratio is then obtained by r(m) := d(m)/(1 − d(m)). Other components of IR follow the scheme of IS, and IR is schematized in Algorithm 8. We point out that IR does not have two variants as IS does, as the definition for E in Example B is not usable because of the unknown normalizing constant for probability p(E|m). See Section A.3 in Appendix for more explanation.
Algorithm 7 Sequential Neural Ratio p1(θ)← p(θ); for r in 1 to R do
repeat sample θi,θ′i ∼ pr(θ); simulate xi ∼ p(x|θi); until n samples are obtained D ← D ∪ {(θi,xi)}ni=1 D′ ← D′ ∪ {(θ′i,xi)}ni=1 train dφ(θ,x) classifying betweenD andD′ with the loss in Eq. 4;
rφ(θ,x)← dφ(θ,x)1−dφ(θ,x) ; pr+1(θ) ∝ rφ(θ,x) · p(θ);
end for return p̂(θ|xo) = pR+1(θ)
Algorithm 8 Iterative Ratio p1(m)← p(m); for r in 1 to R do
repeat sample mi ∼ pr(m); query the oracle: si ← f(mi); until n samples are obtained D ← D ∪ {(mi, si)}ni=1; construct D̃ with m in D satisfying E ; construct D̃′ from p(m); train dφ(m) classifying between D̃ and D̃′; rφ(m)← dφ(m)1−dφ(m) ; pr+1(θ) ∝ rφ(m) · p(m);
end for return {m : (m, s) ∈ D}
Interestingly, a recent SNR-like work, EG-LF-MCMC (Begy & Schikuta, 2021), proposes to train the classifier on tuples of (θ, = ‖x− xo‖) instead of (θ,x). This algorithm can be also seen as a more precise analogy of our Iterative Ratio in the LFI context, as we have pointed out in Section 2.2 that a conceptual link could be drawn between s and −‖x− xo‖.

heading: 3.4 COMPOSITE PROBABILISTIC METHODS
text: Building on the above analogies and framework, we move beyond the above algorithms in this subsection. The previous lines of approach – the direct posterior modeling methods in Section 3.1 and the indirect methods in Section 3.2 and 3.3 – both have their own advantages and disadvantages. The former methods may fail to get accurate inference result due to unmatched proposals, while the latter ones would need extra large amount of computation for the MCMC sampling process before obtaining accurate posterior samples, etc. Here we study composite algorithms that combine the aforementioned ingredients through the lens of our proposed unified framework. Our goal is to combine the strengths from both kinds of methods.
We first introduce Iterative Posterior Scoring (IPS) method illustrated in Algorithm 9. IPS also uses a neural network f̂φ(m) to model the forward mechanism as IS does, and again we use q̃ here to denote the target posterior p(m|E). Instead of applying computational MCMC steps here, we train a second parametrized model qψ to model q̃ by minimizing the KL divergence between them. Notice that the optimization of qψ is restricted within a neural network parameterization family. As a result, in the next round, we can directly utilize qψ(m) to serve as a flexible generative proposal of pr+1(m). Like IS, the IPS algorithm also has two different variants with regard to different choices of E , we name them to be IPS-A and IPS-B. Two choices differ in the detailed construction of distribution q̃(m).
In a similar spirit, we propose the Iterative Posterior Ratio (IPR) algorithm (see Algorithm 10). IPR is close to the IR algorithm in many aspects, but also adopts a second neural network model qψ(m) like IPS. IPR works similarly to IR, in that we also construct q̃(m), taking advantage of rφ(m) and the prior p(m) simply via q̃(m)← rφ(m) · p(m). Note that the usage of two models in IPR is not exactly the same as in IPS: qψ(m) is also achieved via minimizing KL divergence with q̃(m), but the training of model dφ(m) is closer to that in IR rather than the f̂φ(m) in IS. Another similarity between IPR and IR is that IPR also only has one variant, since the Example B is not applicable for this ratio modeling approach (see Section A.4 in Appendix for more details).
Algorithm 9 Iterative Posterior Scoring p1(m)← p(m); for r in 1 to R do
repeat sample mi ∼ pr(m); query the oracle: si ← f(mi); until n samples are obtained D ← D ∪ {(mi, si)}ni=1; fit f̂φ(m) with D; construct q̃(m) with f̂φ(·) and p(m); qψ ← arg minqDKL(q̃(m)‖q); pr+1(m)← qψ(m);
end for return {m : (m, s) ∈ D}
Algorithm 10 Iterative Posterior Ratio p1(m)← p(m); for r in 1 to R do
repeat sample mi ∼ pr(m); query the oracle: si ← f(mi); until n samples are obtained D ← D ∪ {(mi, si)}ni=1; construct D̃ with m in D satisfying E ; construct D̃′ from p(m); train dφ(m) classifying between D̃ and D̃′; rφ(m)← dφ(m)1−dφ(m) ; construct q̃(m) with rφ(m) and p(m); qψ ← arg minqDKL(q̃(m)‖q); pr+1(m)← qψ(m);
end for return {m : (m, s) ∈ D}

heading: 4 EXPERIMENTS

heading: 4.1 SETUP
text: In this section, we systematically evaluate the proposed methods and baselines on four different in-silico biological sequence design benchmarks. In every round, we allow each algorithm to query the black-box oracle for a batch of n sequences mi to obtain their true scores si, with n = 100 for all experiments. The total number of rounds differs across different tasks.
We experiment with our six proposed methods: Iterative Scoring (-A/B) from Section 3.2, Iterative Ratio from Section 3.3, Iterative Posterior Scoring (-A/B) and Iterative Posterior Ratio from Sec-
tion 3.4. Apart from these proposed methods, we consider as baselines a battery of existing methods designed for batched black-box sequence design tasks: (1) Random, a method that randomly select proposal sequences at every round; (2) FB-VAE (Gupta & Zou, 2019) depicted in Section 3.1; (3) Evolution based (Brindle, 1980; Real et al., 2019) sequence design algorithm; and (4) DbAS (Brookes & Listgarten, 2018), described in Section 3.1.
We evaluate these sequence design algorithms by the average score of the top-10 and top-100 sequences in the resulting datasetD at each round. We plot the average score curves with regard to the number of rounds. We also use the area under the curve as a scalar metric for sample efficiency to compare the methods being evaluated. Specifically, since the area depends on the choice of x-axis, we simply cumulate the scores of all rounds to calculate the area. The result could be non-positive, since the score can take negative values.

heading: 4.2 RESULTS
text: Transcription factor binding sites (TfBind). Protein sequences that bind with DNA sequences to adjust their activity are called transcription factors. In Barrera et al. (2016), the authors measure the binding properties between a battery of transcription factors and all possible length-8 DNA sequences through biological experiments. Concretely, we choose 15 transcription factors to serve as 15 different tasks. For each transcription factor, the algorithm needs to search for sequences that maximize the corresponding binding activity score. The size of the search space is |V|L = 48 = 65536. The number of total rounds is fixed to 10. For validation, we follow Angermüller et al. (2020b) and use one task (ZNF200 S265Y R1) for hyperparameter selection. Then we test the algorithms’ performance on the other 14 held-out tasks.
Figure 1 displays a comparison for all ten methods on three of the chosen binding affinity tasks. We can observe that after 10 rounds, our proposed methods perform consistently better than the
backward modeling methods like DbAS and FB-VAE in terms of both Top-10 and Top-100 scores. Among all baselines, the evolution method is the strongest one, and it beats IR on some of the tasks (see complete results in Table 4 and 5 in Appendix). We also find that IS-A and IS-B increase top scores slightly faster than other methods, especially on PBX4 REF R2 and CRX E80A R1. This indicates that composite methods’ way of using parameterized models to replace computational procedures is not the optimal solution for small-scale tasks.
5’ untranslated regions (UTR). The translation efficiency is mainly determined by the sequence of 5’ UTR (Alipanahi et al., 2015). In Sample et al. (2019), the authors create a library of gene sequences with ribosome loading level as labels. They further train a convolutional neural network with this library to predict the relationship between a 5’UTR sequence and the corresponding gene expression level. We use this neural network as an oracle for this benchmark. The length of the gene sequences is fixed to 50, and thus the size of the search space is 450. For this 5’UTR benchmark, we also allow each algorithm to explore for 10 rounds. Figure 2 (left) shows that our proposed composite methods significantly outperform other methods on the UTR task. Different from the results on the TfBind task, forward modeling methods do not achieve the best performance.
Antimicrobial peptides (AMP). Protein modeling has recently become a popular sub-area of machine learning research. We are tasked to generate AMP sequences, which are short protein sequences against multi-resistant pathogens. We train a binary classifier model to classify whether a short protein sequence belongs to AMP and defer the related details to Appendix. This is the only task we consider regarding sequence design with alterable lengths, where the length of sequences ranges from 12 to 60. Since each entry of protein sequence has |V| = 20 different choices on amino acids, the size of search space is ∑60 L=12 20
L. We set the number of total rounds to be 15 for this task. Figure 2 (middle) clearly shows that the performances of IPS-A and IPS-B dominate the AMP generation task, which demonstrates the effectiveness of our composite strategy.
Fluorescence proteins (Fluo). As another protein engineering task, we consider the optimization over fluorescent proteins, which is a commonly used test bed of modern molecular biology. This task is similar to the AMP task introduced above, but its ground-truth measurement relies on a regressor. We use a pretrained model taken from Rao et al. (2019) to act as our task oracle, which is a regressor trained to predict log-fluorescence intensity value over approximately 52,000 protein sequences of length 238 (Sarkisyan et al., 2016). More concretely, the regressor is fit on a small neighborhood of parent green fluorescent protein, and is then evaluated on a more distant protein. The training data is derived from the naturally occurring GFP in Aequorea victoria. The task Fluo’s search space is 20238 and we set the number of rounds to 20. We demonstrate the Fluo results in Figure 2 (right), where we can see both composite methods and DbAS achieve much better results than other approaches. This might signify that for long sequence design tasks, backward modeling is superior to other modeling methods, which is not consistent with LFI (Papamakarios et al., 2019). We also summarize the sample efficiency results for the latter 3 benchmarks in Table 3 and 6, from which we can see that our proposed three composite methods perform remarkably promising results.

heading: 5 CONCLUSION
text: We propose a probabilistic framework that unifies likelihood-free inference and black-box optimization for designing biological sequences. This unified perspective enables us to design a variety of novel composite probabilistic sequence design methods combining the best of both worlds. Extensive experiments have demonstrated the benefits of the unified perspective. While the composite probabilistic methods usually outperform other baseline methods in most sequence design tasks we consider in this work, the key contribution of our paper is not just about the superiority of those composite methods, as different specific tasks might prefer different algorithmic configurations due to no free lunch theorem (Wolpert & Macready, 1997). Actually, we would like to attribute the strong performance to the unified probabilistic framework, which enables us to develop a richer algorithm pool, based on which we can design performant algorithms for particular sequence design tasks.

heading: ACKNOWLEDGEMENT
text: The authors would like to thank Christof Angermueller, Yanzhi Chen, Michael Gutmann, and anonymous reviewers for helpful feedbacks. Jie Fu thanks Microsoft Research Montreal for funding his postdoctoral position at University of Montreal and Mila. Yoshua Bengio acknowledges the funding from CIFAR, Samsung, IBM and Microsoft. Aaron Courville thanks the support of Samsung, Hitachi and CIFAR.

heading: A MORE ABOUT METHODOLOGY
text: A.1 BACKWARD MODELING OF THE MECHANISM
About SMC-ABC. In Algorithm 1, when we pick the top-m θ at the r-th round, the picked parameters actually follow such a distribution
p(θ|xo) ∝ r∑ l=1 pl(θ) · p(‖x− xo‖ < | θ) = ( r∑ l=1 pl(θ) ) · p(‖x− xo‖ < | θ) (5)
where the here is implicitly defined by how ”top” the selection process is, namely the ratio mnr . As a result, we point out that in Algorithm 1, another more complicated form of the last step pr+1(θ)← qφ(θ) is pr+1(θ) ∝ qφ(θ)p(θ)/ ∑r l pl(θ), where an additional renormalizing term is involved. We ignore this term and used the simpler alternative in order to keep our main text clean. We refer interested readers to Beaumont et al. (2009) for more details.
About Sequential Neural Posterior. Define p(x) = ∫ p(x|θ)p(θ)dθ and p̃(x) = ∫ p(x|θ)p̃(θ)dθ for any arbitrary proposal distribution p̃(θ) which is not necessary to be the prior p(θ). What’s more, we define p(θ|x) = p(x|θ)p(θ)/p(x) and p̃(θ|x) = p(x|θ)p̃(θ)/p̃(x) to be the true posterior and the proposal posterior.
Starting from the goal of approximating the true posterior,
arg min q Ep(x)[DKL(p(θ|x)‖q(θ|x))] = arg max q
∫ p(x)dx ∫ p(θ|x) log q(θ|x)dθ
= arg max q
∫ p(θ,x) log q(θ|x)dθx
= arg max q
Ep(θ,x)[log q(θ|x)].
It seems that we can directly train the parameterized neural density estimator qφ in a data driven manner via maxφ ∑ i log qφ(θi|xi) where i is the index for data sample. When the number of training samples as well as the parameterization family of φ are large enough, the obtained qφ would be close enough to the true posterior. This would require the data samples to follow (θi,xi) ∼ p(θ,x) = p(θ)p(x|θ). However, practically one uses a proposal p̃(θ) to first generate some {θi}i and then generate {xi}i by simulation. When the proposal distribution p̃(θ) is not exactly the prior distribution p(θ), the resulting qφ(·|·) would be:
p̃(θ|x) = p(θ|x) p̃(θ)p(x) p(θ)p̃(x) ∝ p(θ|x) p̃(θ) p(θ) , (6)
which is a biased estimation and is not what we want.
Three variants of SNP take different approaches to try to fix this bias. SNP-A (Papamakarios & Murray, 2016) first fits the biased proposal posterior p̃(θ|x) in the aforementioned way and utilize the relation in Eq. 6 to solve for an unbiased estimation. This approach is restricted to mixture of Gaussian distribution family and thus has limited expressiveness. SNP-B (Lueckmann et al., 2017) uses importance sampling to address this issue via maxφ E(x,θ)∼p(x|θ)p̃(θ) [ p(θ) p̃(θ) log qφ(θ | x) ] . One
downside of this approach is the high variance involved by the importance weights p(θ)/p̃(θ). SNP-C (Greenberg et al., 2019) proposes to use reparameterize the proposal posterior by setting q̃φ(θ|x) = qφ(θ|x) p̃(θ)p(θ) 1 Zφ(x) where Zφ(x) = ∫ qφ(θ|x) p̃(θ)p(θ)dθ is the corresponding normalizing factor for x. SNP-C then maximizes E(x,θ)∼p(x|θ)p̃(θ) [log q̃φ(θ | x)]. About Design by Adaptive Sampling. We aim to approximate the posterior via minimizing the KL divergence:
arg min q DKL(p(m|E)‖q(m)) = arg max q
∫ p(m|E) log q(m)dm
= arg max q
∫ p(E|m)p(m) log q(m)dm
= arg max q
Eq̃(m) [ p(m)
q̃(m) p(E|m) log q(m)
] ,
where q̃(m) could be any distribution of m. Brookes et al. (2019) takes this formulation. Brookes & Listgarten (2018) only differs in the place that it ignores the denominator term. According to Angermüller et al. (2020b), we choose the latter variant as one of our baselines because it is more stable in practice. We refer interested readers to Brookes & Listgarten (2018) for more details.
Notice that we are not doing exactly the same things for LFI and black-box sequence design. Since LFI models flexible posterior p(θ|x) which is a distribution for arbitrary x, we can also choose to model p(m|s) for arbitrary s. Nevertheless, in the neural network modeling, conditioning by a scalar value is not an effective approach as the effect of low dimensional scalar value conditioning may be covered by other high dimensional input. Therefore, we choose to directly model the target posterior p(m|E) with a single neural network.
A.2 FORWARD MODELING OF THE MECHANISM
We still use p̃(θ) to denote an arbitrary proposal distribution and p̃(θ,x) := p(x|θ)p̃(θ). Then we have
arg min q Ep̃(θ) [DKL (p(x|θ)‖q(x|θ))] = arg max q
∫ p̃(θ)dθ ∫ p(x|θ) log q(x|θ)dx
= arg max q
Ep̃(θ,x) [log q(x|θ)] .
We point out that with much enough data and large enough expressiveness of the neural density estimator parameterization family, no matter what proposal p̃(θ) is used to provide training samples {(θi,xi)}i ∼ p̃(θ,x), we have the resulting qφ̂(θ|x) equals true likelihood p(x|θ) in the support of the proposal. What SNL gives is an unbiased estimation and thus does not have the same problem as SNP.
Now we elaborate the construction of q̃(m) in Iterative Scoring algorithm. Here we use the notation q̃(m) to denote our approximation of the posterior p(m|E). Notice that we want q̃(m) ∝ p(m) · p(E|m). If we choose Example A to serve as the definition of event E , then the samples of q̃(m) can be obtained in this way: (1) sample m from prior p(m) and (2) accept this sample if f̂φ(m) is larger than threshold s, or otherwise reject it. Alternatively, if we choose Example B, we have q̃(m) ∝ p(m) · exp(f̂φ(m)/τ). Similar to SNL, we do MCMC sampling from this unnormalized probability function.
A.3 MODELING A PROBABILITY RATIO
About Sequential Neural Ratio. Dataset D is generated in the way that (1) first sample θ ∼ p(θ) and (2) simulate x ∼ p(x|θ). Consequently, D follows the distribution p(θ)p(x|θ) = p(θ,x). On the other hand, the other dataset D′ generates θ and x in parallel and independent manner. Notice here θ ∼ p(θ) and x follows the marginal distribution: x ∼ p(x) = ∫ p(θ)p(x|θ)dθ.
Proof of Proposition 1.
Proof. We define a functional F to be the optimization objective:
F [d] = Ea∼p0(a)[log d(a)] + Ea∼p1(a)[log(1− d(a))]
We calculate its functional derivative. For arbitrary function u and infinite small
F [d+ u]−F [d] = Ep0 [log(1 + u
d )] + Ep1 [log(1 + −u 1− d )]
= ∫ u · ( p0 d + −p1 1− d ) +O( )
⇒ lim →0
F [d+ u]−F [d]
= ∫ u · ( p0 d + −p1 1− d ) = ∫ u · δF .
We set the functional derivative to zero:
δF = 0⇒ p0 p1 = d 1− d
⇒ d(a) = p0(a) p0(a) + p1(a) .
This optimal function d∗ apparently takes value in [0, 1].
About Iterative Ratio. Notice that in Algorithm 8 we construct two datasets: D̃ and D̃′. To generate D̃, we need to be able to pick some sequence samples m from D and make the selected ones follow the posterior p(m|E). This procedure will depend on our choice of event E . For Example A this is easy, since we just need to filter out the sequences whose oracle value is smaller than the threshold. However, for Example B, it is hard to do similar things, since given score value from D we only know the unnormalized value of posterior probability, and cannot determine which sequence should be filtered out. The construction of D̃′ which follows prior distribution p(m) is trivial.
A.4 COMPOSITE PROBABILISTIC METHODS
For IPS, the optimization with regard to the second parameterized model qψ(m) is
qψ = arg min q DKL(q̃(m)‖q) = arg max q
∫ q̃(m) log q(m)dm
= arg max q
∫ p(m|E) log q(m)dm = arg max
q
∫ p(E|m)p(m) log q(m)dm.
The exact value of p(E|m) depends on different choices of configuration feature E in Section 2.2. On the other hand, IPR, like IR, also only has one variant, which is with Example A:
qψ = arg min q DKL(q̃(m)‖q) = arg max q
∫ rφ(m)p(m) log q(m)dm,
which is a tractable optimization problem. Both IPR and IR are not fit for Example B since it cannot provide an exact probability value and thus cannot be adopted to construct D̃.

heading: B MORE ABOUT EXPERIMENTS
text: Random method uses no neural network model. FB-VAE uses a VAE model. The encoder of the VAE first linearly transform one-hot input into a hidden feature which is 64 dimension, and then separately linearly transform to a 64-dimension mean output and 64-dimension variance output. The decoder contains a 64×64 linear layer and a linear layer that maps the hidden feature to categorical output. All other methods utilize bi-directional long short-term memory model (BiLSTM) (Hochreiter & Schmidhuber, 1997) with a linear embedding layer. Both the embedding dimension and the hidden size of LSTM is set to 32. For composite methods that use two models, we use one-layer LSTM for each of them. For the other algorithms that only use one LSTM, we set its number of
layers to be two. No Dropout (Srivastava et al., 2014) is used in LSTM models. In this way, the number of parameters of the VAE is slightly larger than that of the two layer BiLSTM, and all methods (except Random) share similar model parameter size.
All the experiments are repeated with fifty random seeds and report the mean value (and also standard deviation in the figure plots). We set p(m) to be uniform prior for all tasks for simplicity, which uniformly samples from the dictionary V for each entry of the sequence. For length alterable task, we first uniformly sample the length between minimum length and maximum length and then sample each entry.
We explain details about evolution based method mentioned in the main text, which can be seen as a substantial example of directed evolution (Chen & Arnold, 1991). Like other model based methods, Evolution also trains an LSTM regressor to predict the score of a sequence, which is further used to assist in the reproduce procedure. The Evolution algorithm maintains a generation list through the whole exploration process. In each round, the method mutates and reproduces the sequences to enlarge the generation list, and then utilizes the learned regressor to select top sequences for the next generation.
For validation, we follow (Angermüller et al., 2020b) and sweep each algorithm for fifty trials and pick the best configuration. We tune learning rate and whether to re-initialize the optimizer for each new round for all methods. We tune threshold for DbAS, FB-VAE and the methods that is with Example A. For the other choice of E , we tune the temperature. For evolution, we tune the number of offsprings for each sequence in generation list, the probability of substitution, insertion and deletion. For TfBind we use ZNF200 S265Y R1 for validation. For UTR, AMP and Fluo, since we only have one oracle instance for each benchmark, we do not use a hold-out validation method.
For TfBind benchmark, we use the following transcription factor instances and treat them as different black-box optimization tasks: ZNF200 S265Y R1, POU6F2 REF R1 8, KLF11 R402Q R1, EGR2 R359W R1, HOXD13 S316C R1, HOXB7 K191R R1, PBX4 REF R2, GFI1B A204T R1, FOXC1 REF R1, KLF1 REF R1, SIX6 REF R1, ARX L343Q R2, CRX E80A R1, ESX1 K193R R1 and VSX1 G160D R1. We do not do post-processing such as score normalization whitening for the data for simplicity. We first calculate the area under curve to summarize the performance in a scalar output, and put the ranking result for each algorithm in Table 4 and Table 5, which provide more details for Table 2. To further investigate the diversity of different algorithms, we choose two TfBind instances (POU6F2 REF R1 and KLF11 R402Q R1) and visualize the resulting sequences with T-SNE (van der Maaten & Hinton, 2008) in Figure 3. We provide two visualization views for both task instances: (1) we uniformly sample 20 sequences from the whole n · R sequences for each algorithm and visualize them; (2) for each algorithm, we visualize 20 sequences uniformly sampled from the last batch (i.e., at the last round). This is notated with “final sequences” in the figure. We do not visualize all the sequences for simplicity. We use Hamming distance in the computation of T-SNE. From Figure 3, we can see that there is no obvious difference for the evaluated methods. This indicates that our proposed methods can achieve better performance while maintaining on-par diversity level with the baselines. This is not exactly consistent to the findings of Angermüller et al. (2020a), which claims some algorithms such as DbAS achieve very limited diversity. We do not use the “optima fraction” metric in Angermüller et al. (2020a;b), since this metric may not deal with multimode oracle landscape well and needs extra unstable computation such as clustering. Besides, this metric cannot generalize to other benchmarks.
We elaborate the construction of our AMP oracle. We use the AMP dataset from (Witten & Witten, 2019) which contains 6,760 AMP sequences. A multilayer perceptron classifier is trained to predict if a protein sequence can prohibit the growth of a particular pathogen in that AMP dataset. This classifier operates on the features extracted by ProtAlbert (Elnaggar et al., 2020) model. Following the setup in (Angermüller et al., 2020b), we treat the predicted logits as the ground-truth measurement. Moreover, we demonstrate the area under Top-10 curves for UTR, AMP and Fluo benchmarks in Table 6, which is a good complement for Table 3 but is missing due to limited space in the main text.

heading: C RELATED WORKS AND DISCUSSION
text: Likelihood-free inference. We have already introduced the main classes of likelihood-free inference algorithms in the main text: (1) Approximate Bayesian Computation (ABC) method (Beaumont et al., 2009; Blum, 2009; Marin et al., 2012; Lintusaari et al., 2017) in Section 3.1; (2) Posterior modeling method that is also stated in Section 3.1, including classical ones (Tran et al., 2015; Li et al., 2017; Chen & Gutmann, 2019) and modern SNP methods (Papamakarios & Murray, 2016; Lueckmann et al., 2017; Greenberg et al., 2019); (3) Likelihood modeling method described in Section 3.2, also containing various classical algorithms (Wood, 2010; Mengersen et al., 2012; Drovandi et al., 2018) and modern SNL variants (Lueckmann et al., 2018; Papamakarios et al., 2019); and (4) Probability ratio modeling methods mentioned in Section 3.3 diverge in estimating likelihood ratio (Gutmann & Hyvärinen, 2010; Gutmann et al., 2018; Brehmer et al., 2020) or likelihood-to-evidence ratio (Thomas et al., 2016; Izbicki et al., 2014), where the latter paradigm is a good fit for LFI problem (Hermans et al., 2019). Besides, there are also works about how to construct low-dimensional summary statistics for LFI (Fearnhead & Prangle, 2012; Chan et al., 2018; Chen et al., 2021).
Machine learning based drug design. Generative modeling and discriminative modeling are two basic ways of thinking in machine learning. In literature for sequence design, generative modeling is also known as cross entropy method. This is a famous kind of design method that is close to our “backward modeling of the mechanism” approach. Cross entropy methods seek to solve an expectation maximization problem (i.e., maxp Ep(m)[f(m)]) where the sequences follow a distribution p. This can also be related to simulated annealing, a large family of black-box optimization algorithm – the sequential neural posterior could be thought to maintain a distribution which is gradually becoming sharper to a delta distribution at the optimal value. On the other hand, we think of this as a way for modeling the posterior p(m|E) and develop corresponding analysis under the probabilistic framework, which is like a more accurate version of cross entropy method. Many related methods (including the ones stated in Section 3.1) train the distribution by likelihood maximization for sequences with large scores, or use some sort of reweighting to achieve similar effects (Rubinstein & Kroese, 2004; de Boer et al., 2005; Neil et al., 2018; Gupta & Zou, 2019; Brookes et al., 2019).
Discriminative modeling usually goes in a “model-based optimization” way (terminology from Angermüller et al. (2020a)), i.e., use a discriminative model f̂(m) to fit the real oracle f(m) and act as a surrogate for it. The surrogate model can replace the true oracle f(m) which involves costly biological experiments. This corresponds to our “forward modeling of the mechanism” in Section 3.2. Bayesian optimization (Shahriari et al., 2016) is a classical example, which utilizes f̂ (typically a Gaussian process model) to define an acquisition function to guide the exploration and exploitation. Many modern biochemical methods also belong to this category (Gómez-Bombarelli et al., 2018; Hashimoto et al., 2018; Yang et al., 2019; Wu et al., 2019; Sample et al., 2019; Liu et al., 2020).
There seems not much related literature about probability ratio estimation based method in this topic. On the other hand, Hashimoto et al. (2018) shares a classification based approach with IR but they differ on how to use the classifier. This algorithm utilizes the learned classifier to update the proposal with multiplicative weights algorithm, making the proposal to have large probability where the classifier logit is small. Other categories of drug design methods include evolution algorithms (Brindle, 1980; Wierstra et al., 2008; Salimans et al., 2017; Yoshikawa et al., 2018; Jensen, 2019; Real et al., 2019; Ahn et al., 2020) that search over the target space with genetic operators like insert, mutation, and crossover, and reinforcement learning (Guimaraes et al., 2017; Neil et al., 2018; Zhou et al., 2019; Shi et al., 2020; Angermüller et al., 2020b) which see the formation of a drug as a Markov decision process and train the policy to learn highly-rewarding drugs. We do not find other work that is similar to our probability ratio modeling approach (Section 3.3) from the literature, which we take as a novel contribution.

Review summaries:
The paper draws on connections between likelihood-free inference and black-box optimization to propose new black-box optimization methods. In general, the goal here is to not find the exact optimum of the black-box objective, but to sample from a set of sequences with high-quality objective. This is akin to the problem of collecting posterior samples in a likelihood-free inference problem.   The paper provides a number of proposed methods and compares them on some benchmark sequence optimization problems that have appeared in recent literature.

The authors describe a mapping from likelihood free inference to black-box sequence optimization, then use this mapping to link common algorithms in both fields. They go on to describe novel black-box sequence design algorithms induced by known LFI algorithms. Empirical results show their methods are competitive on standard datasets.

In this paper the authors draw direct parallels between likelihood-free inference (LFI) and black-box sequence design.  This allows that authors to draw parallels between existing methods from the LFI and black-box sequence design literatures.  In a few cases there is no direct analog in the black-box sequence design literature for a given LFI algorithm, and so the authors are able to immediately propose such an algorithm.  The authors also present a number of "composite" methods that combine ideas from a number of these approaches.