Generating Focussed Molecule Libraries for Drug Discovery with Recurrent Neural Networks

1.7 MB·

by  Marwin H. S. Segler

Segleretal.
RESEARCH
Generating Focussed Molecule Libraries for Drug
Discovery with Recurrent Neural Networks
Marwin HS Segler1*, Thierry Kogej2, Christian Tyrchan3 and Mark P Waller4
7
1
0 Abstract
2 In de novo drug design,computationalstrategies are used to generate novel molecules with good affinityto the

n desired biological target. In this work, we show that recurrent neural networks can be trained as generative
a
models for molecular structures, similar to statistical language models in natural language processing. We
J
demonstrate that the properties of the generated molecules correlate very well with the properties of the
5
molecules used to train the model. In order to enrich libraries with molecules active towards a given biological

] target, we propose to fine-tune the model with small sets of molecules, which are known to be active against
E
that target.
N
Against Staphylococcus aureus, the model reproduced 14% of 6051 hold-out test molecules that medicinal
. chemists designed, whereas against Plasmodium falciparum (Malaria) it reproduced 28% of 1240 test
s
c molecules. When coupled with a scoring function, our model can perform the complete de novo drug design
[ cycle to generate large sets of novel molecules for drug discovery.

1 Keywords: computer-assisted drug design; recurrent neural networks
v
9
2
3 1 Introduction lions of existing or billions of virtual molecules.[5]
1
Chemistry is the language of nature. Chemists speak Searching can be carried out using similarity-based
0
it fluently and have made their discipline one of the metrics, which provides a quantifiable numerical indi-
.
1 true contributors to human well-being, which has cator of closeness between molecules. In contrast, in
0
“change[d] the way you live and die”.[1] This is par- de-novo drug design,one aims to directly create novel
7
1 ticularly true for medicinal chemistry. However,creat- molecules that are active towards the desired biologi-
: ing novel drugs is an extraordinarily hard and com- caltarget.[6,7]Here,likeinanymoleculardesigntask,
v
plex problem.[2] One of the many challenges in drug the computer has to
i
X design is the sheer size of the search space for novel i create molecules,
r molecules. It has been estimated that 1060 drug-like ii score and filter them, and
a moleculescouldpossiblybesyntheticallyaccessible.[3]
iii searchforbettermolecules,buildingontheknowl-
Chemists have to select and examine molecules from
edge gained in the previous steps.
this large space to find molecules that are active to-
Task i, the generation of novel molecules, is usu-
wards a biological target. Active means for example
ally solved with one of two different protocols.[7]
that a molecule binds to a biomolecule, which causes
One strategy is to build molecules from predefined
an effect in the living organism,or inhibits replication
groupsofatomsorfragments.Unfortunately,theseap-
of bacteria. Modern high-throughput screening tech-
proachesoften leadto molecules that arevery hardto
niques allow to test molecules in the order of 106 in
synthesise.[8] Therefore,anotherestablished approach
the lab.[4] However, larger experiments will get pro-
is to conduct virtual chemical reactions based on ex-
hibitively expensive. Given this practical limitation of
pert coded rules, with the hope that these reactions
in vitro experiments, it is desirable to have compu-
could then also be applied in practice to make the
tational tools to narrow down the enormous search
molecules in the laboratory.[9]These systems giverea-
space. Virtual screening is a commonly used strat-
sonabledrug-likemolecules,andareconsideredas“the
egy to search for promising molecules amongst mil-
solution” to the structure generation problem.[2] We
*Correspondence:[e-mail] generally share this view. However, we have recently
1InstituteofOrganicChemistry&CenterforMultiscaleTheoryand
shown that the predicted reactions from these rule-
Computation,WestfälischeWilhelms-Universität,Münster,Germany
Fulllistofauthorinformationisavailableattheendofthearticle based expert systems can sometimes fail.[10] Also, fo-
Segleretal. Page2of17
cussing ona smallset ofrobustreactionscanunneces- similar to the training data. Intuitively, a generative
sarily restrict the possibly accessible chemical space. model for molecules trained on drug molecules would
Task ii, scoring molecules and filtering out unde- "know" how valid and reasonable drug-like molecules
sired structures, can be solved with substructure fil- look like, and could be used to generate more drug-
ters for undesirable reactive groups in conjunction like molecules. However, for molecules, these models
with established approaches such as docking[11] or have been studied rarely, and rigorously only with
machine-learning(ML) approaches.[7, 12,13] The ML traditional models such as Gaussian mixture models
approaches are split into two branches: Target predic- (GMM).[37, 41] Recently, recurrent neural networks
tion classifies molecules into active and inactive, and (RNNs) have emerged as powerful generative models
quantitative structure-activity relationships (QSAR) in very different domains, such as natural language
seek to quantitatively predict a real-valued measure processing,[42] speech,[43] images,[44] video,[45] for-
for the effectiveness of a substance (as a regression mal languages,[46] computer code generation,[47] and
problem). As molecular descriptors, Signature Fin- musicscores.[48]Inthiswork,wehighlighttheanalogy
gerprints, Extended-Connectivity (ECFP) and atom of language and chemistry, and show that RNNs can
pair (APFP) fingerprints and their fuzzy variants are also generate reasonable molecules. Furthermore, we
thede-facto standardtoday.[14–16]ConvolutionalNet- demonstratethatRNNscanalsotransfertheirlearned
works on Graphs are a recent addition to the field knowledgefromlargemoleculesetstodirectlyproduce
of molecular descriptors.[17, 18] Random Forests and novel molecules that are biologically active by retrain-
Neural Networks are currently the most widely used ing the models on small sets of already known actives.
machine learning models for target prediction.[19–31] We test our models by reproducing hold-out test sets
This leads to task iii, the search for molecules with of known biologically active molecules.
therightbindingaffinitycombinedwithoptimalmolec-
2 Methods
ular properties. In earlier work, this was performed
(amongothers)with classicalglobaloptimisationtech- 2.1 Representing Molecules
niques, for example genetic algorithms or ant-colony To connect chemistry with language, it is important
optimisation.[7, 32] Furthermore,de novo designis re- to understand how molecules are represented. Usu-
lated to inverse QSAR.[33–36] While in de novo de- ally,theyaremodeledbymoleculargraphs,alsocalled
sign design, a regular QSAR mapping X → y from Lewis structures in chemistry. In molecular graphs,
molecular descriptor space X to properties y is used atoms are labeled nodes. The edges are the bonds be-
as the scoring function for the global optimizer, in tween atoms, which are labeled with the bond order
inverse QSAR one aims to find an explicit inverse (e.g. single, double or triple). One could therefore en-
mapping y → X, and then maps back from optimal vision having a model that reads and outputs graphs.
points in descriptor space X to valid molecules. How- Several common chemistry formats store molecules in
ever, this is not well defined, because molecules are such a manner. However, in models for natural lan-
inherently discrete. Several protocols have been devel- guage processing, the input and output of the model
opedtoaddressthis,forexampleenumeratingallstruc- areusuallysequencesofsingleletters,stringsorwords.
tures within the constraintsof hyper-rectanglesin the We therefore employ the Smiles format, which en-
descriptor space.[33–38] Gómez-Bombarelli et al. pro- codes molecular graphs compactly as human-readable
posedtolearncontinuousrepresentationsofmolecules strings. Smiles is a formal grammar which describes
with variational auto-encoders, based on the model molecules with an alphabet of characters,for example
by Bowman et al.,[39] and to perform Bayesian op- candCforaromaticandaliphaticcarbonatoms,Ofor
timisation in this vector space to optimise molecular oxygen, -, = and # for single, double and triple bonds
properties.[40]Nevertheless,this approachwasnotap- (see Figure1).[49]To indicaterings,anumber is intro-
plied to create active drug molecules, and did not suc- duced at the two atoms where the ring is closed. For
ceedinoptimisingmorecomplexmolecularproperties, example, benzene in aromatic Smiles notation would
such as emission color and delayed fluorescence decay be c1ccccc1.Side chains are denoted by roundbrack-
rate (kTADF).[40] ets. To generate valid Smiles, the generative model
In this work, we suggest a novel, completely data- would have to learn the Smiles grammar, which in-
drivende novo drug design approach.It relies only on cludeskeepingtrackofringsandbracketstoeventually
a generative model for molecular structures, based on close them. In morphine, a complex natural product,
a recurrent neural network, that is trained on large the number of steps between the first 1 and the sec-
sets of molecules. Generative models learn a proba- ond 1, indicating a ring, is 32. Having established a
bility distribution over the training examples; sam- link between molecules and (formal) language,we can
pling from this distribution generates new examples now discuss language models.
Segleretal. Page3of17
O x1:n = (x1,...,xn) and an initial state vector h0, and
N N O returns a sequence of state vectors h1:n = (h1,...,hn)
OH
O N N and a sequence of output vectors y1:n = (y1,...,yn).
The RNN consists of a recursively defined function R,
Caffeine Ibuprofen
CN1c2ncn(C)c2C(=O)N(C)C1=O CC(C)Cc1ccc(cc1)C(C)C(O)=O which takes a state vector hi and input vector xi+1
HO andreturnsanew statevectorhi+1.Anotherfunction
O maps a state vector hi to an output vector yi.[50]
O H
N
HO
Morphine RNN(h0,x1:n)=h1:n,y1:n (2)
[H][C@]12C=C[C@H](O)[C@@H]3Oc4c5c(C[C@H]1N(C)CC[C@@]235)ccc4O
hi =R(hi−1,xi) (3)
Figure1 Examplesofmoleculesandtheir Smiles
representation. Tocorrectly createsmiles,themodelhasto yi =O(hi) (4)
learnlongtermdependencies, forexampletocloserings
(indicated bynumbers) andbrackets. The state vector hi stores a representationofthe in-
formationaboutallsymbolsseeninthesequencesofar.
As analternativetothe recursivedefinition,the recur-
2.2 Language Models and Recurrent Neural Networks
Given a sequence of words (w1,...,wi), language mod- a) xt b) x1 x2 x3 x4
els predict the distribution of the (i + 1)th word
wi+1.[50]Forexample,ifalanguagemodelreceivesthe ht-1 R,O ht h0 R,O h1 R,O h2 R,O h3 R,O h4
θ θ θ θ θ
sequence "Chemistry is", it would assign different
probabilities to possible next words. "fascinating", yt y1 y2 y3 y4
"important", or "challenging" would receive high Figure 2 a)RecursivelydefinedRNNb)ThesameRNN,
probabilities,while "runs"or"potato"wouldreceive unrolled. Theparameters θ (the weightmatricesoftheneural
network) aresharedoveralltimesteps.
verylowprobabilities.Languagemodels canbothcap-
ture the grammatical correctness ("runs" in this sen-
tence is wrong) and the meaning ("potato" does not rent network can also be unrolled for finite sequences
make sense). Language models are implemented for (seeFigure2).AnunrolledRNNcanbe seenasavery
example in message autocorrection in many modern deep neural network, in which the parameters θ are
smartphones. Interestingly, language models do not shared among the layers, and the hidden state ht is
have to use words. They can also be based on charac- passed as an additional input to the next layer.Train-
ters or letters.[50] In that case, when receiving the se- ing the unrolledRNNto fitthe parametersθ canthen
quence of characterschemistr,it would assigna high simply be done by using backpropagationto compute
probability to y, but a low probability to q. To model the gradients with respect to the loss function, which
molecules instead of language, we simply swap words is categoricalcross-entropyin this work.[50]
or letters with atoms, or, more practically, characters As the specific RNN function, in this work, we use
in the Smiles alphabet, which form a (formal) lan- the Long Short Term Memory (LSTM), which was in-
guage.Forexample,ifthemodelreceivesthesequence troduced by Hochreiter and Schmidhuber.[51] It has
c1ccccc,thereisahighprobabilitythatthenextsym- been used successfully in many natural language pro-
bol would be a "1", which closes the ring, and yields cessing tasks,[42] for example in Google’s Neural Ma-
benzene.
chine Translation system.[52] For excellent in-depth
Moreformally,toasequenceS ofsymbolssi atsteps discussions of the LSTM, we refer to the articles by
ti ∈T, the language model assigns a probability of Goldberg,[50]Graves,[53]Olah,[54]andGreffet al.[55]
T To encode the Smiles symbols as input vectors xt
, we employ the "one-hot" representation.[53] This
Pθ(S)=Pθ(s1)·YPθ(st|st−1,...,s1) (1) means if there are K symbols, and k is the symbol
t=2
to be input at step t, then we can construct an input
where the parameters θ are learned from the training vector xt with length K, whose entries are all zero ex-
set.[50] In this work, we use a recurrent neural net- cept the k-th entry, which is one. If we assume a very
work (RNN) to estimate the probabilities of Equation restrictedsetofsymbols{c,1,\n},inputc wouldcor-
1. In contrast to regular feedforward neural networks, respond to xt = (1,0,0), 1 to xt = (0,1,0) and \n to
RNNs maintain state, which is needed to keep track xt =(0,0,1).
of the symbols seen earlier in the sequence. In ab- TheprobabilitydistributionPθ(st+1|st,...,s1)ofthe
stractterms,anRNNtakesasequenceofinputvectors next symbol given the already seen sequence is thus
Segleretal. Page4of17
Step t 1 2 3 ... 7 8
s c 1 c ... c 1
t
xt (1,0,0) (0,1,0) (1,0,0) (1,0,0) (0,1,0)
h h h h h
1 2 3 8 9
c1ccccc1\n
yt (0.2,0.8,0) (0.9,0,0.1) (0.8,0.1,0.1) (0.2,0.8,0) (0.1,0.1,0.9) gather
symbols
s 1 c c ... 1 \n
t+1
Figure3 TheSymbolGenerationandSamplingProcess.Westartwitharandomseedsymbols1,herec,whichgetsconverted into
aone-hotvectorx1 andinputintothemodel.Themodelthenupdates itsinternal stateh0 toh1 andoutputs y1,whichisthe
probabilitydistributionoverthenextsymbols.Here,samplingyieldss2=1.Convertings2 tox2,andfeedingittothemodelleadsto
updated hidden stateh2 andoutput y2,fromwhichcansampleagain.Thisiterativesymbol-by-symbolprocedure canbecontinued
aslongasdesired. Inthisexample,westopitafterobservinganEOL(\n)symbol,andobtaintheSmilesforbenzene. Thehidden
statehi allowsthemodel tokeeptrackofopened brackets andrings,toensurethatthey willbeclosedagainlater.
a multinomial distribution, which is estimated using 2.3 Transfer Learning
the output vector yt of the recurrent neural network For many machine learning tasks, only small datasets
at time step t by areavailable,whichmightleadtooverfittingwithpow-
erfulmodelssuchasneuralnetworks.Inthissituation,
exp(yk) transfer learning can help.[58] Here, a model is first
t
Pθ(st+1 =k|st,...,s1)= K exp(yk′) (5) trained on a large dataset for a different task. Then,
Pk′=1 t the modelisretrainedonthe smallerdataset,whichis
also called fine-tuning. The aim of transfer learning is
where yk corresponds to the k-th element of vector to learngeneralfeatures on the biggerdata set, which
t
yt.[53] Sampling from this distribution would then also might be useful for the secondtask in the smaller
allow generating novel molecules: After sampling a data regime. To generate focussed molecule libraries,
Smiles symbol st+1 for the next time step t+1, we we first train on a large, general set of molecules,
can construct a new input vector xt+1, which is fed then perform fine-tuning on a smaller set of specific
into the model, and via yt+1 and Equation 5 yields molecules,andafterthatstartthesamplingprocedure.
Pθ(st+2|st+1,...,s1). Sampling from the latter gener-
ates st+2, which serves again also as the model’s in- 2.4 Target Prediction
put for the next step (see Figure 4). This symbol-by- To verify whether the generated molecules are active
symbol sampling procedure is repeated until the de- on the desired targets,standard target prediction was
sired number of characters has been generated.[53] employed. Machine learning-based target prediction
aims to learn a classifier c : M → {1,0} to decide
To indicate that a molecule is "completed", each
whetheramoleculem∈moleculardescriptorspaceM
molecule in our training data finishes with an “end
is activeornotagainstatarget.[12,13]The molecules
of line” (EOL) symbol, in our case the single charac-
are split into actives and inactives using a threshold
ter \n (which means the training data is just a simple
Smiles file). Thus, when the system outputs an EOL, on a measure for the substance effectiveness. pIC50 =
−log10(IC50) is one of the most widely used metrics
a generated molecule is finished. However, we simply
continue sampling, thus generating a regular Smiles for this purpose. IC50 is the half maximal inhibitory
concentration,thatistheconcentrationofdrugthatis
file that contains one molecule per line.
requiredtoinhibit50%ofabiologicaltarget’sfunction
In this work, we used a network with three stacked in vitro.
LSTM layers, using the keras library.[56] The model To predict whether the generated molecules are ac-
was trained with back propagation through time,[53] tive towards the biological target of interest, target
using the ADAM optimizer at standard settings.[57] prediction models (TPMs) were trained for all the
Tomitigatetheproblemofexplodinggradientsduring tested targets (5-HT2A, Plasmodium falciparum and
training, a gradient norm clipping of 5 is applied.[53] Staphylococcus aureus).We evaluated Random Forest,
Segleretal. Page5of17
LogisticRegression,(Deep)NeuralNetworksandGra- where n = |GN ∩ T| is the number of reproduced
dient Boosting Trees (GBT) as models with ECFP4 molecules from T by sampling a set GN of |GN| = N
(Extended Connectivity Fingerprint with a diameter molecules from the fine-tuned generative model, and
of 4) as the molecular descriptor.[14, 15] We found m=|RM ∩T| is the number of reproduced molecules
thatGBTsslightlyoutperformedallothermodels,and from T by sampling a set RM of |RM|=M molecules
used these as our virtual assay in all studies (see Sup- from the generic, unbiased generative model trained
porting Information for details). ECFP4 fingerprints only on the large dataset. Intuitively, EOR indicates
weregeneratedwithCDKversion1.5.13.[59,60]Scikit- how much better the fine-tuned models work when
Learn,[61] xgBoost[62] and keras[56] were used as the compared to the general model.
machine learning libraries. For 5-HT2A and Plasmod-
ium, molecules are considered as active for the TPM 3 Results and Discussion
if their IC50 reportedin ChEMBL is < 100nm, which In this work, we address two points: First, we want
translates to a pIC50 > 7, whereas for Staphylococcus, to generate large sets of diverse molecules for vir-
we used pMIC > 3.
tual screening campaigns. Second, we want to gener-
ate smaller, focussed libraries enriched with possibly
2.5 Data
active molecules for a specific target. For the first
ThechemicallanguagemodelwastrainedonaSmiles
task, we can train a model on a large, general set
filecontaining1.4millionmoleculesfromtheChEMBL of molecules to learn the Smiles grammar. Sampling
database, which contains molecules and measured
fromthismodelwouldgeneratesetsofdiverse,butun-
biological activity data. The Smiles strings of the
focused molecules. To address the second task, and to
molecules were canonicalized (which means finding
obtain novel active drug molecules for a target of in-
a unique representation that is the same for isomor-
terest, we performtransferlearning:We selecta small
phic molecular graphs)[63, 64] before training with
set of known actives for that target and we refit our
the CDK chemoinformatics library, yielding a Smiles
pre-trained chemical language model with this small
file thatcontainedonemoleculeperline.[59,60]Ithas
data-set. After eachepoch, we sample from the model
to be noted that ChEMBL contains many peptides,
to generate novelactives. Furthermore, we investigate
natural products with complex scaffolds, Michael ac-
if the model actually benefits from transfer learning,
ceptors, benzoquinones, hydroxylamines, hydrazines
by comparing it to a model trained from scratch on
etc. which is reflected in the generatedstructures (see
the small sets without pre-training.
below).Thiscorrespondsto72millionindividualchar-
acters, with a vocabulary size of 51 unique characters.
3.1 Training the recurrent network
51 characters is only a subset of all Smiles symbols,
We employed a recurrent neural network with three
since the molecules in ChEMBL do not contain many
stacked LSTM layers, each with 1024 dimensions,
oftheheavyelements.Aswehavetosetthenumberof
and each one followed by a dropout[65] layer, with
symbols as a hyperparameter during model construc-
a dropout ratio of 0.2, to regularise the neural net-
tion,andthemodelcanonlylearnthedistributionover
work. The model was trained until convergence,using
the symbols present in the training data, this implies
abatchsizeof128.TheRNNwasunrolledfor64steps.
thatonlymoleculeswiththese51Smilessymbolsseen
It had 21.3×106 parameters.
during training can be generated during sampling.
During training, we sampled a few molecules from
The 5-HT2A, the Plasmodium falciparum and the
themodelevery1000mini-batchestoinspectprogress.
Staphylococcus aureus datasetwerealsoobtainedfrom
Within a few 1000 steps, the model starts to output
ChEMBL. The molecules for the hold-out test sets
valid molecules (see Table 1).
were removed from the training data.
3.2 Generating Novel Molecules
2.6 Model Evaluation
To generate novel molecules, 50,000,000 Smiles sym-
To evaluate the models for a test set T, and a set of
bols were sampled from the model symbol-by-symbol.
molecules GN generated from the model by sampling,
we report the ratio of reproduced molecules |GN∩T|, This correspondedto 976,327 lines, from which 97.7%
|T| were valid molecules after parsing with the CDK
and enrichment over random (EOR), which is defined
toolkit. Removing all molecules already seen during
as,
training yielded 864,880 structures. After filtering out
duplicates,weobtained847,955novelmolecules.Afew
n generated molecules were randomly selected and de-
EOR= |GN| (6) picted in Figure 4. The Supporting Information con-
m
|RM| tains more structures. The created structures are not
Segleretal. Page6of17
Table 1 Moleculessampledduringtraining.
Batch Generated Example valid
0 Oc.BK5i%ur+7oAFc7L3T=F8B5e=n)CS6RCTAR((OVCp1CApb) no
1000 OF=CCC2OCCCC)C2)C1CNC2CCCCCCCCCCCCCCCCCCCCCCC no
2000 O=C(N)C(=O)N(c1occc1OC)c2ccccc2OC yes
3000 O=C1C=2N(c3cc(ccc3OC2CCC1)CCCc4cn(c5c(Cl)cccc54)C)C yes
S O O
OONO NHOHN O HN NF NH2N OH NO S NON HNO OO 1105 t-SNE of Phys.Chem Descriptors CGheEneMrBatLed
N N HNN NNO FFF Cl ON FNNN AdNHONH OHN ON O 05
O Cl
Figure4 Afewrandomlyselected, generatedmolecules.Ad= −5
Adamantyl
−10
−15
−20 −15 −10 −5 0 5 10 15 20
just formallyvalid,but arealsomostly chemicallyrea-
sonable.
Figure 5 t-SNEprojection of7physicochemical descriptors of
In order to check if the de novo compounds could randommoleculesfromChEMBL(blue)andmolecules
be considered as valid starting points for a drug dis- generated withtheneuralnetwork trainedonChEMBL
(green), totwounitlessdimensions. Thedistributions ofboth
covery program, we applied the internal AstraZeneca
setsoverlapsignificantly.
filters.[66]AtAstraZeneca,thisflaggingsystemisused
to determine if a compound is suitable to be part of
the high-throughputscreening collection (if flaggedas
“core” or “backup”) or should be restricted for particu- Furthermore, we analysed the Bemis-Murcko scaf-
lar use (flagged as “undesirable” since it contains one folds of the training molecules and the sampled
or several unwanted substructures, e.g. undesired re- molecules.[70]Bemis-Murckoscaffoldscontainthering
active functional groups). The filters were applied to systems ofa molecule and the moieties that link these
the generated set of 848 k molecules and they flagged ring systems, while removing any side chains. They
most of them, 640 k (75%), are either core or backup. represent the scaffold, or “core” of a molecule, which
Since the same ratio (75%) of core and backup com- series of drug molecules often have in common. The
pounds hasbeenobservedforthe ChEMBLcollection, number of common scaffolds in both sets, divided by
wethereforeconcludethatthealgorithmgeneratespre- the union of all scaffolds in both sets (Jaccard index)
ponderantly valid screening molecules and faithfully is 0.12, which indicates that the language model does
reproduces the distribution of the training data. not just modify side chainsubstituents, butalso intro-
To determine whether the properties of the gener- duces modifications at the molecular core.
ated molecules match the properties of the training
data from ChEMBL, we followed the procedure of 3.3 Generating Active Drug Molecules and Focused
Kolb:[67] We computed several molecular properties, Libraries
namely molecular weight, BertzCT, the number of H- 3.3.1 Targeting the 5-HT2A receptor
donors, H-acceptors, and rotatable bonds, logP and To generate novel ligands for the 5-HT2A receptor,
total polar surface area for randomly selected sub- we first selected all molecules with pIC50 > 7 which
sets from both sets with the RDKit[68] library ver- weretestedon5-HT2A fromChEMBL(732molecules),
sion 2016.03.1.Then, we performed dimensionality re- andthenfine-tunedourpre-trainedchemicallanguage
duction to 2D with t-SNE (t-Distributed Stochastic model on this set. After each epoch, we sampled
NeighborEmbedding,atechniqueanalogousto PCA), 100,000 chars, canonicalised the molecules, and re-
which is shown in Figure 5.[69] Both sets overlap al- moved any sampled molecules that were already con-
most completely, which indicates that the generated tained in the training set. Following this, we evalu-
moleculesverywellrecreatethepropertiesofthetrain- ated the generated molecules of each round of retrain-
ing molecules. ing with our 5-HT2A target prediction model (TPM).
Segleretal. Page7of17
In Figure 4, the ratio of molecules predicted to be 3.3.2 Targeting Plasmodium falciparum (Malaria)
active by the TPM after each round of fine-tuning Plasmodium falciparum is a parasite that causes the
is shown. Before fine-tuning (corresponding to epoch most dangerous form of Malaria.[71] To probe our
0), the model generates almost exclusively inactive model on this important target, we used a more
molecules. Already after 4 epochs of fine-tuning the challenging validation strategy. We wanted to inves-
model produced a set in which 50% of the molecules tigate whether the model could also propose the same
are predicted to be active. moleculesthatmedicinalchemistschosetoevaluatein
publishedstudies.Totestthis,first,theknownactives
againstPlasmodiumfalciparumwithapIC50 >8were
selected from ChEMBL. Then, this set was split ran-
1.0 Influence of fine-tuning epochs on the ratio of actives domly into a training (1239 molecules) and a test set
(1240 molecules). The chemical language model was
then fine-tuned on the training set. 7500 molecules
0.8
were sampled after each of the 20 epochs of refitting.
Ratio of actives00..46 TEOa#bRle:E2pnrIRCice5hp0mroednuTtcrtoaivninegr.krnanoTwdoensmta.ctivGeesni.nmthoels.PlasRmeopdroiudm. teEstOsRet.
1 >8 1239 1240 128,256 28% 66.9
0.2 2 >8 100 1240 93,721 7% 19.0
3 >9 100 1022 91,034 11% 35.7
0.0
0 5 10 15 20 This yielded128,256unique molecules.Interestingly,
Epochs of fine-tuning
we found that our model was able to "redesign" 28%
Figure6 Epochs offine-tuningvsratioofactives. of the unseen molecules of the test set. In compari-
sontomoleculessampledfromtheunspecific,untuned
model, an Enrichment over Random (EOR) of 66.9 is
obtained.Withasmallertrainingsetof100molecules,
Diversity Analysis In order to assess the novelty of the model can still reproduce 7% of the test set, with
the de novo molecules generated with the fine-tuned an EOR of 19.0. To test the reliance on pIC50 we
model,anearestneighborsimilarity/diversityanalysis chose to use another cut-off of pIC50 > 9, and took
has been conducted using a commonly used 2D fin- 100 molecules in the training set and 1022 in the test
gerprint(ECFP4) basedsimilarity method (Tanimoto set. 11% of the test set could be recreated, with an
index).[67]Figure7showsthedistributionofthenear- EOR of 35.7. To visually explore how the model pop-
est neighbor Tanimoto index generated by comparing ulates chemical space, Figure 8 shows a t-SNE plot of
all the novel molecules and the training molecules be- the ECFP4 fingerprintsofthe testmoleculesand2000
generated molecules that were predicted to be active
fore and after n epochs of fine-tuning. For each bin,
by the target prediction model for Plasmodium falci-
the white bars indicate the molecules generated from
parum.Itindicatesthatthemodelhasgeneratedmany
the unbiased, general model, while the darker bars in-
similar molecules around the test examples.
dicatethemoleculesafterseveralepochsoffine-tuning.
Within the bins corresponding to lower similarity, the
3.3.3 Targeting Staphylococcus aureus (Golden Staph)
numberofmoleculesdecreases,whilethebinsofhigher
To evaluate a different target, we furthermore con-
similarity get populated with increasing numbers of
ducted a series of experiments to reproduce knownac-
molecules. The plot thus shows that the model starts
tivemoleculesagainstStaphylococcus aureus.Here,we
to output more and more similar molecules to the
used actives with a pMIC > 3. MIC is the Mean In-
target-specifictrainingset.Notably,afterafewrounds
hibitory Concentration, the lowest concentration of a
of training not only highly similar molecules are pro-
compound that prevents visible growth of a microor-
duced, but also molecules covering the whole range of ganism. As above, the actives were split into a train-
similarity, indicating that our method could not only ing and a test set. However, here, the availability of
deliver close analogs but new chemotypes or scaffold the data allows larger test sets to be used. After fine-
ideas to a drug discovery project.[5] To have the best tuning on the training set of 1000 molecules (Table 3,
of both worlds,that is diverse and focussedmolecules, Entry1),ourmodelcouldretrieve14%ofthe6051test
we therefore suggest to sample after each epoch of re- molecules.Whenscalingdowntoasmallertrainingset
training and not just after the final epoch. of 50 molecules (the model gets trained on less than
Segleretal. Page8of17
Figure7 Nearest-neighbour Tanimotosimilaritydistribution ofthegenerated moleculesfor5-HT2A afternepochsoffine-tuning
againsttheknown actives.Thegenerated moleculesaredistributed overthewholesimilarityrange.Generated moleculeswitha
mediumsimilaritycanbeinteresting forscaffold-hopping.[5]
Table 3 Reproducting known activesintheStaphylococcus testset.EOR:Enrichment overrandom.
Entry pMIC Train. Test Gen.mols. Reprod. EOR
1 >3 1000 6051 51,052 14% 155.9
2 >3 50 7001 70,891 2.5% 21.6
3a >3 50 7001 85,755 1.8% 6.3
4b >3 50 7001 285 0% —
5c >3 0 7001 60,988 6% 59.6
aFine-tuning learningrate=10−4.bNoPretraining. c8Generate-Test cycles.
1%ofthedata!),itcanstillreproduce2.5%ofthetest point for the typical cyclical drug discovery process:
set, and performs 21.6 times better than the unbiased Molecules get designed, synthesised, and then tested
model (Table 3, Entry 2). Using a lower learning rate in assays. Then, the best molecules are selected, and
(0.0001,Entry3)forfine-tuning,whichisoftendonein based on the gained knowledge new molecules are de-
transferlearning,doesnotworkaswellasthestandard signed,whichclosesthecycle.Therefore,asafinalchal-
learning rate (0.001, Entry 2). We additionally exam- lenge for our model, we simulated this cycle by iterat-
inedwhetherthemodelbenefitsfromtransferlearning. ing molecule generation("synthesis"), selection of the
When trained from scratch,the model performs much best molecules with the machine learning-based tar-
worsethanthepretrainedandsubsequentlyfine-tuned getprediction("virtualassay")andretrainingthelan-
model (see Figure 9 and Table 3, Entry 4). Pretrain- guage model with the best molecules ("design") with
ingonthelargedatasetisthuscrucialtoachievegood Staphylococcus aureus as the target. We thus do not
performance against Staphylococcus aureus. use asetofknownactivestostartthe structuregener-
ation procedure (see Figure 10).
3.4 Simulating Design-Synthesis-Test Cycles Westartedwith100,000sampledmoleculesfromthe
The experiments we conducted so far are applicable if unbiased chemical language model. Then, using our
one already knows several actives. However, in drug target prediction model, we extracted the molecules
discovery, one often does not have such a set to start classified as actives. After that, the RNN was fine-
with. Therefore, high throughput screenings are con- tuned for 5 epochs on the actives, sampling ≈10,000
ducted toidentify afew hits,whichserveasastarting molecules after each epoch. The resulting molecules
Segleretal. Page9of17
t-SNE of ECFP4
20 GPleansemroadteiudm, P fr.e Tdeicstte Sde Atctives 1) RRNetNrain Target Prediction 1) Score
15 2) Sample Model 2) Filter Inactives
10
5
New Active Molecules
0
−5 Figure 10 Schemeofourdenovo designcycle.Moleculesare
generated bythechemicallanguagemodel andthenscored
−10 withthetargetprediction model (TPM).Theinactivesare
filteredout, andtheRNNisretrained. Here,theTPMisa
−15
machinelearningmodel,butitcouldalsobearobot
−20 conducting synthesisandbiological assays,oradocking
−20 −15 −10 −5 0 5 10 15 20
program.
Figure8 t-SNEplotofthepIC50>9testset(blue)andthede
novo moleculespredicted tobeactive(green). Thelanguage
model populateschemical spacearound thetestmolecules. 3.5 Why does the model work?
Our results presented in Section 3.2 show that the
general model trained on a large molecule set has
learned the Smiles rules and can output valid, drug-
Fine-Tuning vs Direct Training like molecules,whichresemblethe trainingdata.How-
1.4 frtoympe scratch ever, sampling from this model does not help much if
pretrained
we want to generate actives for a specific target: We
1.2
would have to generate very large sets to find actives
1.0 for that target among the diverse range of molecules
Test Loss0.8 tshcoermesoidnelocurreeaxtepse,rwimheicnhtsi.sindicatedbythehighEOR
0.6
0.4
0.2 0.16
0.0 0.14
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Epochs
0.12
FaFuiingreeu-urtseund9iantDgaistffeheterwepnirttehttrr1aa0iinn0e0indgtmrasotinrdaientlgegpaieenrsfdoor6nm0s5th1beettSetstetarpethxhyaalmoncptolrecasci.nuisng Frequency00..0180
fromscratch(lowertestloss[crossentropy] isbetter). med
Nor0.06
0.04
0.02
werefilteredwiththetargetpredictionmodel,andthe
0.00
newactivesappendedtothe activesfromthe previous 0 10 20 30 40 50 60
Levenshtein Distance
round, closing the loop.
Figure 11 HistogramofLevenshtein (Stringedit)distances of
Alreadyafter8iterations,themodelreproduced416 theSmilesofthereproduced moleculestotheirnearest
of the 7001 test molecules from the previous task, neighbour inthetrainingset(Staphylococcus aureus,model
retrained on50actives).Whileinmanycasesthemodel
whichis6%(Table3,Entry5),andexhibitsandEOR makeschangesofafewsymbolsintheSmiles,resembling
of 59.6. This EOR is higher than if the model is re- thetypicalmodificationsappliedwhenexploringseriesof
compounds, thedistribution ofthedistancesindicatesthatthe
trained directly on a set of 50 actives (Entry 2). Ad-
RNNalsoperformsmorecomplex changesbyintroducing
ditionally, we obtained 60,988 unique molecules that largermoietiesorgenerating moleculesthatarestructurally
different, butisofunctional tothetrainingset.
the target prediction model classified as active. This
demonstrates that in combination with a target pre-
diction or scoring model, our model can also perform When fine-tuned to a set of actives, the probability
the complete de novo-design cycle. distributionoverthemoleculescapturedbyourmodel
Segleretal. Page10of17
is shifted towards molecules active towards our tar- a specific biologicaltarget, which enables the creation
get.Tostudythis,wecomparetheLevenshtein(String of novel molecules with the desired activity. By iter-
edit)distanceofthegeneratedSmilestotheirnearest ating cycles of structure generationwith the language
neighboursinthetrainingsetinFigure11.TheLeven- model, scoring with a target prediction model (TPM)
shtein distance of e.g. benzene c1ccccc1and pyridine and retraining of the model with increasingly larger
c1ccncc1 would be 1. Figure 11 shows that while the sets of highly scored molecules, we showedthat we do
model often seems to have made small replacements notevenneedasetofactiveknownactivemoleculesto
in the underlying Smiles, in many cases it also made start our procedure with, as the TPM could also be a
more complex modifications or even generated com- docking program, or a robot conducting synthesis[72]
pletely differentSmiles.Thisis supportedalsoby the and biological testing.
distribution of the nearest neighbour fingerprint simi-
We see three main advantages of our method. First,
laritiesoftrainingandrediscoveredmolecules(ECFP4,
it is conceptually orthogonal to established molecule
Tanimoto, Figure 12). Many rediscovered molecules
generationapproaches,as it learnsa generativemodel
are in the medium similarity regime.
for molecular structures. Second, our method is very
Because we perform transfer learning, during fine-
simple to setup, train and to use, and can be adapted
tuning, the model does not "forget" what it has
to different datasets without any modifications to the
learned.Aplausibleexplanationwhy themodelworks
model architecture, and does not depend on hand-
is therefore that it can transfer the modifications that
encoded expert knowledge. Furthermore, it merges
areregularlyappliedwhenseriesofmoleculesarestud-
structuregenerationandoptimisationinonemodel.A
ied, to the molecules it has seen during fine-tuning.
weakness of our model is interpretability. In contrast,
existing de-novo design methods settled on virtual re-
actionstogeneratemolecules,whichhasadvantagesas
1.0 Nearest Neighbour Distribution - Training vs Reproduced Molecules it minimises the chance of obtaining "overfit", weird
molecules, and increases the chances to find synthesiz-
0.8 able compounds.[2, 7]
To extend our work, it is just a small step to cast
0.6 molecule generation as a reinforcement learning prob-
milarity lem, where the pre-trained LSTM generator could be
Si seen as a policy, which can be encouraged to create
0.4
better molecules with a rewardsignalobtainedfroma
target prediction model.[73] In addition, different ap-
0.2
proaches for target prediction, for example docking,
could be evaluated.[7, 11]
0.0
DeepLearningisnotapanacea,andwejoinGawehn
et al. in expressing “some healthy skepticism” regard-
Figure12 Violinplotofthenearest-neighbour ing its application in drug discovery.[26] Generating
ECFP4-Tanimoto similaritydistribution ofthe50training
molecules that are almost right is not enough, be-
moleculesagainsttherediscovered moleculesinTable3,Entry
2.Thedistribution suggests themodel haslearnedtomake cause in Chemistry, a miss is as good as a mile, and
typicalsmallfunctional groupreplacements, butcanalso drug discovery is a “needle in the haystack” problem
reproduce moleculeswhicharenottoosimilartothetraining
data. – in which also the needle looks like hay. Neverthe-
less, given that we have shown in this work that our
model can rediscover those needles, and other recent
developments,[26, 74–76] we believe that deep neu-
4 Conclusion
ral networks can be complimentary to established ap-
Inthis work,wehaveshownthatrecurrentneuralnet-
proachesindrugdiscovery.Thecomplexityoftheprob-
worksbasedontheLongShortTermMemory(LSTM)
lem certainly warrants the investigation of novel ap-
canbe appliedto learna statisticalchemicallanguage
proaches. Eventually, success in the wet lab will de-
model. The model can generate large sets of novel
termine if the new wave[21] of neural networks will
molecules with similar physico-chemical properties to
prevail.
the training molecules. This can be used to generate
librariesforvirtualscreening.Furthermore,wedemon-
stratedthatthemodelperformstransferlearningwhen
Competinginterests
fine-tuned to smaller sets of molecules active towards Theauthorsdeclarethattheyhavenocompetinginterests.