Latent Molecular Optimization for Targeted Therapeutic Design

The goal of computational drug design is to assist in the discovery of molecules with therapeutic potential. Given the vast and rapidly expanding knowledge from molecular biology, one may know the cellular machinery (usually a protein) responsible for a given disease, and can thus devise a targeted therapeutic, designed to affect a particular biomolecule within the target pathway. This paradigm, known as molecular targeting therapy, has shown promise in therapeutic development (e.g. cancer [1, 2]). However, discovering and constructing a molecule (or ligand) capable of binding a given target is non-trivial. The space of small molecules with therapeutic potential is estimated to be on the order of $10^{60}$ [3, 4]. More problematic is the discreteness of chemical structure space, in which minor perturbations can lead to large changes in biochemical properties (such as affinity), preventing continuous optimization approaches. Further, strong binding affinity to the target is merely necessary, not sufficient, for a useful therapeutic: the presence of off-target effects, cost, and pharmacological ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties also play a role in the efficacy or viability of a drug [5, 6]. Therapeutic design is therefore a multi-objective optimization problem of considerable difficulty, causing the pharmaceutical industry to struggle to find viable ligands able to reach the end of the development pipeline [7, 8].

One promising route to alleviating this problem is the use of computational methods, such as virtual screening [9], molecular docking [10], and quantitative structure-activity relationship (QSAR) models [11, 12], all of which attempt to predict or characterize affinity of a small molecule to a target. More recently, the advent of deep learning has improved various prediction tasks in cheminformatics [13, 14], such as bioactivity scoring [15, 16], organic synthesis planning [17], and toxicity estimation [18]. Of particular interest, however, is the de novo design of molecules with desirable properties [19, 20] (also known as the inverse QSAR problem [21, 22, 23]), which has been approached in a variety of ways, including Bayesian optimization in latent space [24, 25, 26, 27, 28], reinforcement learning [29], deep generative models [30, 31], and recurrent neural networks (RNNs) [32, 33]. Many of these methods, however, do not address the problem of targeted design, while those that do suffer from short-comings, such as relying on specialized biochemical models, requiring bound complexes for affinity estimation, or not utilizing the structural data contained in other protein-ligand pairs.

In this work, we present a novel method to perform targeted design of small molecule agents via optimization in the latent space of a generative model. First, using a dataset of protein-ligand complexes (PLCs), we extract and separate the protein binding sites as graph-valued data and the ligands as discrete chemical SMILES strings [34]. A graph convolutional approach is used to map the protein sites to a vector signature $P$, while the Junction Tree Variational Autoencoder (JTVAE) [26] is used to map the discrete chemical to a latent vector $C$. Several differentiable models are then learned via neural networks: (1) a direct mapper, which attempts to reconstruct $C$ from $P$ alone, (2) a binding affinity estimator, and (3) intrinsic property regressors, which predict the ease of synthesis, drug-likeness, and toxicity of $C$. This pipeline is illustrated in Fig. 1. All these models are utilized in a gradient-based optimization algorithm over latent chemical space, which is given only an embedded target protein site $P$ and is tasked with finding a small molecule with both strong binding affinity to the target as well as desirable pharmacological properties.

Our main contributions are to provide ways to circumvent some of the limitations of current approaches. First, by working in a latent vector space, we can use continuous optimization techniques to perform efficient gradient-based search over chemical structures. Second, by utilizing a generic representation of protein surface sites, we can take advantage of the information held in datasets of bound protein-ligand pairs, meaning we do not necessarily require large amounts of biochemical data for a target protein (i.e. we attempt to generalize properties of other complexes to our target). Third, we avoid docking by learning a differentiable scoring function that takes in a protein surface site and ligand separately. To our knowledge, this is the first work utilizing these techniques for targeted molecular design, which we hope will prove useful for future endeavors in this area.

Given the importance of de novo chemical design to the pharmaceutical development process, a considerable body of work has been devoted to computational methods for improving it (e.g. [14, 20]). Much recent research on constructing continuous embeddings of chemicals has used Bayesian optimization to generate molecules with desirable properties [24, 25, 26, 27, 32, 28]. However, these methods do not perform targeted chemical design, aimed at a particular molecule of interest.

Inspired by the success of stochastic gradient-based optimization in training deep neural networks, we also take advantage of the continuous nature of the latent chemical space, using ADAM [35] for optimization, rather than a Bayesian approach. This has the disadvantage of requiring differentiable models for all of the quantities we wish to optimize; however, deep learning models can be used to satisfy this requirement.

Separately, methods of predicting affinity from structural information have also been recently developed [16, 15]. Such networks take a PLC as input, and use the structural information from the complex to score the affinity of the bound pair. The main difference between our method and these approaches is that our affinity prediction network does not receive structural information on the interaction; i.e. the protein site and putative binding chemical are received independently, as embedded vectors. This makes the estimation problem more difficult, since detailed atomic interaction information is not available to the predictor, but more general, since it can be applied to any protein site-chemical pair. This allows us to avoid a docking step (i.e. to generate a PLC) during our optimization.

The most similar studies to ours perform targeted molecular generation. For instance, the work by Takeda et al. [36] assists de novo design by visualization of various properties over chemical space. Other examples include works using RNNs [29, 33] to generate SMILES strings with desired properties, in particular affinity to a target; these works, as well as the work by Blachke et al. [37], train their models on biological activity data for specific protein targets (i.e. using a specific QSAR model). In other words, these models will work best on target protein molecules that already have considerable data available, concerning known binders and their affinities. This is in contrast to our method, which attempts to perform a generic mapping from any given target protein site (i.e. “patch” of atoms from a protein) to a binding chemical. As such, our approach is to rely on learning a general representation of affinity, based on structural information, rather than specializing to a specific input. Thus, our algorithm can be applied to targets for which no biochemical data is known; the predictions of the model will be based on the patterns learned from other complexes. The price to pay for this generality is the need for structural information; however, the exponential growth of such data bodes well for this paradigm [38]. Of course, when such activity data is available for a target, it should be utilized; such information could easily be integrated into our algorithm, simply by adding another term to the energy function, based on the activity predicted from the specific QSAR model.

Targeted drug design cannot only take predicted binding strength into account; it must also consider molecular properties of a candidate therapeutic. Indeed, molecules with good binding capabilities, but poor pharmacokinetics and toxicity properties, can result in costly failures as drugs [62, 63]. We therefore consider the predicted toxicity and several drug-likeness properties of the chemical as well.

First, we considered estimating the toxicological properties of a given small molecule. We used the Tox21 dataset [64], which includes binary data on the activation of five cellular stress response pathways. Not all chemicals have data for all toxicity measures; hence, we use the average toxicity of the measures that were present as the output label, denoted $L_{\text{tox}}(C)$. However, the dataset is very imbalanced, since most of the molecules have a zero average toxicity (see Supplementary Material). We therefore linearly upweighted positive data points via: $W_{\text{tox}}(C)=(W_{m}-1)L_{\text{tox}}(C)+1$, where we fixed $W_{m}=5$ based on the frequency of positive values in the dataset. For our estimator, let $f_{T}:\mathcal{L}\rightarrow[0,1]$ be the toxicity predictor. Its architecture is the same as that used to compute $p_{B}$, but without batch normalization. Fig. 4 (bottom left) shows the distribution of the error $\text{Err}_{T}=f_{T}(C)-L_{\text{tox}}(C)$ on a held-out test set of size 250 (mean with SE of $|\text{Err}_{T}|$: $0.14\pm 0.009$; PC: $0.54$). Note that the prediction errors are skewed to the positive values, due to the weighting function.

We also looked at three chemical properties, known to be associated with drug-likeness: the partition coefficient logP, which controls the lipophilicity (and thus affects the pharmacokinetic and ADMET properties of the molecule [6, 65]), the Quantitative Estimation of Drug-Likeness (QED), which attempts to quantify the aesthetic judgment of medicinal chemists [66], and the synthetic accessibility score (SAS), which estimates the difficulty in synthesis for a given compound [67]. We use the ZINC250K dataset, as extracted by [24] from the ZINC dataset [68]. Given a molecule $C$, denote the logP, QED, and SAS as the components of a vector $\phi(C)=(\phi_{\text{logP}},\phi_{\text{QED}},\phi_{\text{SAS}})$. We construct a function $f_{P}:\mathcal{L}\rightarrow\real^{3}$ that estimates $\phi(C)$, implemented as a neural network with two hidden layers (sizes 120 and 60), trained using MSE loss, and regularized with dropout ($p=0.5$) and weight decay.

The prediction results are shown as the distribution of relative errors (i.e. $|\phi_{i}(C)-f_{P}(C)_{i}|/\text{mean}(\phi_{i})$) on a held-out test set of size 20000 in Fig. 4 (top row). The mean and SE for logP, SAS, and QED are $0.26\pm 0.002$, $0.1\pm 0.0006$, and $0.14\pm 0.001$ (PC: $0.8$, $0.87$, and $0.37$). This suggests that the estimated values of the various chemical drug-likeness properties are within a reasonable range of their true values.

Given the differentiable models above, it is straightforward to define the gradient descent procedure in the latent space (depicted in Fig. 5). We need only to define the energy we are minimizing:

for a fixed target input $P$ and variable latent chemical $C$, where $f_{P}=(\widehat{\phi}_{\text{logP}},\widehat{\phi}_{\text{QED}},\widehat{\phi}_{\text{SAS}})$ and $f_{B}=(p_{B},\widehat{B}_{\text{DSX}})$. The first term in equation 4 acts as the targeting term, ensuring the resulting chemical has good predicted binding affinity to the target; the second term encourages the molecule to have pharmacologically desirable chemical properties. We also view this latter intrinsic energy as a regularizer, preventing the optimization from exploiting anomalous behaviour of the affinity predictor. We chose $\alpha=(-10,1/200)$ and $\gamma=(-0.5,-1,0.1,1)$.

Two terms have additional transformations applied to them. The first is an offset rectifier on the predicted DSX: $g_{h}(x)=\max(-250,x)$, used to limit the term’s influence when it had started to reach the edge of biophysically plausible values (few real DSX values were below $-300$). The second is a quadratic function on the logP: $g_{q}(x)=-4x(x-5)/25$, which encourages the logP to stay near the middle of $[0,5]$. While sufficiently high logP is required for physiological transport of the drug, a logP $>5$ violates Lipinski’s “rule of five” [69] and a logP $>3$ increases probability of toxicity [62]; as such, some literature suggest aiming for a logP within ${\sim}$1-3 [65].

We tested our optimization algorithm on the target protein sites from the held-out PLC test set of size 1000 utilized above (i.e. they were unknown to any of the models). The resulting components of the energy function over time are shown in Fig. 6. It is clear that the predicted binding probability and DSX are indicating strong estimated affinity for the majority of input targets by the end of the optimization. All of the intrinsic properties are also improved: logP stays near 2, while QED steadily increases, and both toxicity and SAS decrease.

Next, we qualitatively visualize some results. While we cannot say whether the chemicals output by the latent space optimization are good binders in reality or not, we can check if they have chemical moieties in common with the known binding ligand (which neither the algorithm nor the affinity model had access to), which would suggest the model learned a structural pattern in the protein site and associated latent molecule data that generalized to unseen graphs. See Fig. 7 for examples of this phenomenon occurring. For instance, in the first row of the figure, the known binder is close to being a subgraph of the output chemicals, while in the second and third row, the resulting chemical after optimization is qualitatively more similar to the binding ligand.

Finally, we discuss quantitative evaluation of our optimization algorithm (see Fig. 8). This was done by docking generated ligands to their target proteins, and then estimating the DSX score of the bound complex; we compare the results to the DSX scores that result from generating random chemicals and docking one to each protein target. We used Open Babel [70] for preprocessing and rDock [71] for protein-ligand docking. In detail, for each target protein in the test set $P_{i}$, we generate two chemicals: $C_{i}$, the output of our optimization algorithm on $P_{i}$, and $R_{i}$, a ligand drawn from the latent prior of the JTVAE. Using rDock, we then dock each of $R_{i}$ and $C_{i}$ to $P_{i}$ separately, obtaining two new bound complexes, from which we can estimate DSX scores $D_{r}(i)$ and $D_{c}(i)$. For docking, since each $P_{i}$ has a known binder, the “reference ligand” approach to cavity generation is used, which restricts the input chemicals to bind to the same site as the known ligand; each docking is run 5 times, and we take the conformation with minimal DSX. Denote the ordered sets of DSX scores of the true, random, and latent binders as $D_{T}$, $D_{R}$, and $D_{C}$ respectively; we also consider the set of differences $\Delta=\{D_{r}(i)-D_{c}(i)\;\forall\;i\}$, so that a positive difference indicates that the optimized chemical has a better DSX score than its random counterpart. Out of 1000 PLCs in the test set, we obtained 961 pairs of docked complexes. In 631 of these (65.7%), the optimized DSX was lower than its random counterpart; further, while only 246 complexes with random ligands had a DSX less than $-100$, there were 477 optimized ligands with DSX score below that. The DSX differences $\Delta$ had a median of 14.5, mean of 9.3, and standard deviation of 51.5. Altogether, these results suggest the optimized ligands are enriched with better binders than one would expect to obtain by random sampling.

One issue is that, when the docking or scoring algorithm fails, we are essentially comparing two random ligands. If we want to look only at the data for which all pipeline components appear to have produced reasonable results (i.e. to see if it is better than random then), as well as eliminate physically implausible outliers, we can restrict our attention to the case when the DSX is less than $-100$. On this subset, 443 (92.9%) of the optimized docked complexes have better DSX scores than their random counterparts. The differences on this subset ${\Delta}_{S}=\{D_{r}(i)-D_{c}(i)\,|\,D_{c}(i)<-100\}$ had a median of 35.4, mean of 36.4, and standard deviation of 26.3. Overall, while these results rely on simulated dockings, rather than empirical experiment, they still provide independent validation of the utility of our approach in generating targeted chemical ligands.

We have proposed an algorithm for targeted therapeutic design: given an input site on a molecule of interest, we wish to generate an agent with satisfactory predicted binding affinity and chemical properties. While both protein binding sites and ligands are inherently discrete and highly structured entities, the use of a graph convolutional signature extraction technique and a deep generative latent variable encoder over chemical structures allows us to embed them into continuous vector spaces. By learning differentiable models of target affinity and intrinsic molecular properties, we are able to utilize gradient-based optimization methods to perform targeted de novo design.

In terms of future work, one could consider other scoring functions (e.g. based on learning [72]), models (e.g. manifold-theoretic geometric learning techniques [73] or using adversarial generative models), data modalities (e.g. augmenting with biochemical data or QSAR models), and search regularizers (e.g. penalizing latent distance from the origin, to help ensure molecular validity). Additional validation, e.g. via other docking methods or molecular dynamics, would also be useful, particularly given that docking can sometimes be unreliable [74, 75]. Overall, our results suggest that this approach could be a promising avenue for future progress in computational drug design.