In the following study, authors aimed to improve the affinity of the compound (**1**) by substituting the dimethylamino group with a more voluminous 2-(benzylamino) moiety (Giordanetto et al., 2012). The novel compound (**2**) showed improved potency (IC50 = 1.9µM) and efficiency (LE = 0.37 and LLE = 4.52) toward p110β. The rationale for this chemical modification relies on the observation that the bulkier substituents might target amino acid residues M804 and W812 in the proximal selectivity pocket. Finally, compound (**3**) was synthesized by introducing the naphthyl group, which in turn attributed to the nanomolar potency (IC50 = 0.093µM) and improved p110β selectivity profile of compound (**3**).

#### **Identification of mitogen-activated protein kinase-interacting kinase 1 inhibitors**

Mitogen-activated protein (MAP) kinase interacting kinases 1 and 2 (**MNK1** and **MNK2**) carry out phosphorylation reaction of eukaryotic translation initiation factor 4E (eIF4E) on serine 209 (Wendel et al., 2007). This translation factor is involved in different cellular pathways, such as Ras/Raf/MEK/ERK and PI3-kinase/protein kinase B (Akt) signaling pathways (Proud, 2015). The overexpression of phosphorylated eIF4E leads to several malignant diseases, such as lymphomas, breast cancer, and glioblastoma (Astanehe et al., 2012). The significance of MNK1/2 enzymes in malignant transformation of the cell has led to high demand for drug design of MNK1/2 inhibitors.

One remarkable study was performed in 2010, where Oyarzabal et al. identified a highly potent and efficient fragment entirely by in silico modeling. In this comprehensive study, authors combined different virtual screening techniques to identify pharmacological tools for MNK1 inhibition. Initially, the Centro Nacional de Investigaciones Oncológicas (CNIO) database was filtered according to the molecular weight (<300 Da) and calculated solubility values (threshold −4 mol/L). By performing this prefiltering procedure, the authors extracted 42,168 fragment-like compounds for virtual screenings (Oyarzabal et al., 2010).

Availability of the crystal structure of MNK2 complexed with staurosporine (PDB: 2HW7) enabled creating minimal substructure, required for crucial interactions with MNK2 (**Figure 8**). The GOLD software (Jones et al., 1997) used in this study was able to reproduce the binding mode of staurosporine in MNK2. The virtual substructure was docked in the crystal structure of MNK1 (PDB: 2HW6) to similarly elucidate crucial amino acid interactions in the ATP binding pocket. MNK1 pharmacophore prepared in this way was used for pharmacophore fitting study, as a molecular docking alternative and 92 compounds were extracted according to their goodness of fit with the pre-defined substructure.

Structure-based virtual screening protocols were combined with ligand-based virtual screenings of CNIO database and external virtual database of compounds collected by the authors who performed this study. These strategies involved 2D-substructural searches, 2D Tanimoto structural similarity, Feature Trees similarity, and three-dimensional shape and electrostatic similarities based on two reported MNK1 inhibitors. Finally, the authors selected 1,236 compounds for biochemical MNK1 assay and 26 of them were active. The hit ratio of this screening was 2.10% and 10 different scaffolds were represented. Interestingly, one compound (**Figure 8**, compound 29) demonstrated nanomolar MNK1 (IC50 = 646 nM) and MNK2 (IC50 = 575 nM) inhibition. Additionally, at the cellular level, compound 29 showed an antiproliferative effect against acute myeloid leukemia cell line (MV4:11, EC50 = 17µM) with dose-dependent decrease in phosphorylation on serine 209 in eIF4E. In conclusion, this study identified 26 hit molecules as MNK1 inhibitors, with 19 of them as fragments with high ligand efficiency values. Among the 26 identified hits, there were 10 diverse chemotypes represented for further drug design studies.

Researchers from A-STAR were particularly interested in imidazopyridazine scaffold (chemotype III in the study of Oyarzabal et al.) as a starting fragment for lead optimization of MNK1/2 inhibitors (Yang et al., 2018). Extensive SAR study of imidazopyridazine derivatives was based on in silico conclusions defined in their previous computational study (Kannan et al., 2017). Concisely, researchers in this study aimed to modify the

screening are labeled in green; (B) the lead optimization strategies starting from fragment ET-38766 to clinical candidate ETC-206.

heterocyclic core in positions 3 and 6, with later modification of the imidazopyridazine scaffold (swapping cores strategy). All the derivatives synthesized in this study were firstly examined by molecular docking studies in Glide 2017-3 software (www. schrodinger.com). By detailed computational analysis of the important amino acid residues in the ATP pocket of MNK1/2 kinases, researchers performed initial lead optimization of the fragment (compound ET-38766) to compound 27 (**Figure 8B**). Novel compound 27 bears imidazopyrazine scaffold, with improved potency against MNK-1 and MNK-2, cell permeability and improved pharmacokinetic properties. After finding optimal substituents in positions 3 and 6, the final step of lead optimization was focused on detailed DFT study to select the final heterocyclic core of MNK1/2 inhibitors. Initially, it was unclear from molecular dynamics (MD) simulations whether the imidazopyrazine N-7 contributes favorably to the binding affinity of MNK inhibitors. To examine this, the authors performed DFT study and demonstrated that N-7 is mostly solvent exposed, thus the final selected heterocyclic core was imidazopyridine. The most promising compound 48 (**Figure 8B**)

later designated as ETC-206, was presented as superior compared to other derivatives in the study. This compound was investigated for the synergism with dasatinib in vivo and currently is in phase I clinical trial for the blast crisis chronic myeloid leukemia (BC-CML).

#### **Computational approaches in rational discovery of allosteric kinase inhibitors**

Although targeting of highly conserved ATP-binding site by Type I and Type II inhibitors provides limited selectivity, inhibiting multiple kinases with a single small-molecule inhibitor was proven to be a useful strategy for therapeutic intervention. However, development of highly selective small-molecule kinase inhibitors remains a pressing concern where targeting of allosteric sites emerged as a promising approach (Wu et al., 2015a). Some of the advantages of targeting allosteric sites include increased selectivity and low toxicity of such inhibitors due to low evolutional conservation of allosteric sites compared to orthosteric (ATP-binding) sites (Fang et al., 2013). Additionally, overcoming of point mutation-associated drug resistance, especially for mutations in the ATP-binding site reported for almost all of ATP-competitive inhibitors, could be another advantage of developing allosteric kinase inhibitors (Gibbons et al., 2012).

While exploitation of allosteric sites represents a very promising strategy, it remains challenging from the aspect of rational drug discovery. Some of the major obstacles include identification of allosteric binding sites, which are usually hidden in less populated higher energy conformations of the proteins. Those conformations are poorly accessible to current experimental methods of structural biology (Lu et al., 2018a). Additionally, allosteric effectors are susceptible to "mode switching," where minor chemical modification of ligand induces critical change in activity (Wenthur et al., 2014). Although known CADD workflows for discovery of drugs directed to orthosteric binding sites are being used in allosteric inhibitors discovery (Rastelli et al., 2014; Schoepfer et al., 2018), they provide limited utility rising the need for development of more spatialized tools and workflows (Greener and Sternberg, 2018).

Identification of allosteric pockets is a crucial first step in rational discovery of allosteric inhibitors. As will discussed below, a plethora of computationally inexpensive methodologies have been developed for this purpose and many of them are even implemented as web servers. While these methodologies provide fast and inexpensive highway in the discovery of druggable allosteric pockets, proper understanding of the allosteric mechanism is impossible without considering underlying conformational landscape and free-energy profiles where more computationally demanding molecular dynamics based approaches have a predominant role. In this review, we discuss few examples of computational methodologies used for direct discovery of novel allosteric sites and/or allosteric kinase inhibitors. For detailed description of recent breakthroughs in computational methodologies used for allosteric inhibitors discovery in general, the interested reader is referred to the recent reviews (Wagner et al., 2016; Lu et al., 2019).

#### Automatic Computational Tools/Web Servers to Investigate Allostery

Structure-based computational tools AlloSite and recently advanced descendant AllositePro (http://mdl.shsmu.edu.cn/ AST/) are intended for fast detection of allosteric site in input PDB structures. Initial detection of allosteric sites is based on Fpocket, a fast open source protein pocket detection software package based on Voronoi tessellation (Le Guilloux et al., 2009). While Allosite uses a machine-learning model to rerank detected pockets in terms of their allosteric character, AllositePro additionally implements normal-mode analysis (NMA) perturbation with elastic network models to account for protein flexibility. NMA is a technique developed for investigation of the vibrational motion of a harmonic oscillating system in the immediate vicinity of its equilibrium. Under assumption that the potential energy landscape in the vicinity of a minimized atomic structure is approximately harmonic, NMA eliminates the need to integrate the equations of motion and makes NMA much less computationally demanding compared to MD (Bahar and Rader, 2005). Zhang et al. demonstrated utility of AllositePro in identification of novel allosteric site on CDK2 kinase. Existence of novel site was validated in mutagenic analysis (Song et al., 2017). Recently, the same group developed AlloFinder, integrated allosterome mapping, and virtual screening workflow implemented as web server (http:// mdl.shsmu.edu.cn/ALF/). AlloFinder relies on AllositePro algorithm for detection of allosteric sites, Allolike filter for pre-filtering of ligand library to enrich allosteric-like compounds (Wang et al., 2012), AutoDock Vina algorithm for docking (Trott and Olson, 2010), and Alloscore empirical scoring function for scoring allosteric modulator-protein complexes (Li et al., 2016). In the final step, alosterome mapping is used to detect highly similar allosteric sites among known human allosteric sites and to rule out selective ligands. This approach was retrospectively validated on several kinase targets (Huang et al., 2018).

CavityPlus (http://www.pkumdl.cn:8000/cavityplus/index. php) is another web server for detection of potential allosteric sites that works on similar principle (Xu et al., 2018). CavityPlus is aimed to detect potential binding sites on the surface of a given protein and rank them based on ligandability and druggability scores. This server integrates several functionalities: CAVITY for detection and scoring of potential binding sites (Yuan et al., 2013); CavPharmer for generation of receptor-based pharmacophores (Chen et al., 2014); CorrSite for prediction of allostery based on NMA motion correlation analysis between allosteric and orthosteric sites (Ma et al., 2016); CivCys for detection of binding sites for covalent inhibitors (Zhang et al., 2017). Functionalities of CavityPlus were successfully used for identification of allosteric binding site on Polo-like kinase 1 (Plk1). Subsequent molecular-docking-based virtual screening on allosteric site resulted in identification of few potent Plk1 inhibitors (Yun et al., 2016).

Another successful implementation of web server based tools for allosteric drug discovery is Kinase Atlas (https://kinaseatlas.bu.edu/) (Yueh et al., 2019). Kinase Atlas is systematic collection of mostly unexplored allosteric sites (binding hot spots) calculated for 4,910 PDB structures of 376 distinct kinases. The hot spots are identified by FTMap. This method places molecular probes (small organic molecules) on a dense grid around the protein and finds favorable positions using an empirical energy function and CHARMM potential. After clustering of obtained positions for each probe, regions that bind several probe clusters are marked as hot spots (Kozakov et al., 2015). Authors of the study identified novel allosteric site on CDK2 and screened library of 1,280 molecules using disulphide-based fragment screening. Two potent and novel allosteric inhibitors were described.

#### Molecular Dynamics-Based Approaches to Investigate Allostery

Molecular dynamics-based approaches in rational discovery of allosteric kinase inhibitors have potential to provide exclusive insight in atomic-level dynamical mechanism of allostery, to explore conformational landscape and capture kinase conformational states inaccessible to current experimental methodologies. Therefore, molecular-dynamics-based approaches, even though being computationally intensive, could detect previously unknown conformations and hidden allosteric binding pockets (Guo and Zhou, 2016; Lu et al., 2018a).

Combination of conventional MD simulations with other standard SBDD approaches resulted in identification of novel allosteric sites and discovery of novel allosteric ligands in several cases. For example, Perez et al. identified novel inhibitory allosteric site and inhibitors of p38α by using MD simulations starting from the X-ray structure of binary complex of p38α and its interacting partner MAPK-activated protein kinase 2 (MK2). MD simulations permitted definition of pharmacophoric features of small peptide inhibitors derived from sequence of MK2. Subsequent virtual screening study resulted in first small molecule allosteric inhibitor for identified binding site (Gomez-Gutierrez et al., 2016). Cournia et al. verified existence of allosteric site on human PI3Kα previously described in murine PI3Kα using combination of FTMap, MD, and in vitro assays. Intriguingly, MD simulations revealed different binding mode of studied allosteric inhibitor in murine, WT, and mutant forms of PI3Kα and consequent differences in propagation of allosteric signal to orthosteric ATP-binding site (Gkeka et al., 2015).

Computational costs of insufficient conformational sampling often limit application of conventional MD simulations in investigating allostery phenomena. Currently, there is a large gap between the time scale which can be reached in MD simulations and that observed in experiments. Several strategies for enhancing the sampling of MD simulations have been proposed (Aci-Seche et al., 2016; Yang et al., 2019). Two recently reported studies demonstrating full power of enhanced sampling methods (Markov-state modeling based adaptive sampling and parallel tempering in the well-tempered ensemble) are discussed below with special reference to atomic-level description of allosteric communication and discovery of cryptic allosteric pockets.

Pande et al. investigated activation pathway of c-Src kinase using massively distributed MD simulations (550 µs) on Folding@HOME (Shirts and Pande, 2000) Markov-state modeling (MSM) and adaptive sampling algorithms in order to provide description of factors underlying thermodynamics and kinetics of c-Src activation and to identify key structural intermediates (Shukla et al., 2014). Briefly, MSM models represent kinetical description of a system's underlying freeenergy landscape, useful for characterization of probability of dynamical transitions between conformational states identified in many independent MD simulations and for extrapolation of long time system's behavior (Sengupta and Strodel, 2018). In this study intermediate conformational state which could be stabilized to block the c-Src activation pathway, was described through MSM analysis for the first time. Further analysis on identified c-Src conformational state revealed the existence of allosteric pocket and surprisingly high structural similarity to known complex of CDK2 bound to allosteric inhibitor— ANS (Betzi et al., 2011). Further simulations confirmed binding of ANS to the novel allosteric site of c-Src and blockage of activation process by stabilization of intermediate states. Additionally, the long-range residues coupling analysis identified myristate-binding pocket as another potential target for development of allosteric modulators of c-Src. Taken together, results of this study highlighted large-scale MD coupled with MSM modeling as an indispensable tool for identification of novel conformational states, potential allosteric pockets, and study of mechanisms of allostery in kinases.

In another example, authors explored the possibility of bidirectional communication between allosteric so-called PIFpocket and ATP-binding site in PDK1 protein kinase using a combination of experimental techniques and enhancedsampling simulations [parallel tempering simulations in the well-tempered ensemble (PT-WTE)] (Schulze et al., 2016). Results of PT-WTE MD revealed bidirectional mechanisms of communication between the ATP-binding site and allosteric site. Interestingly, this study for the first time demonstrated how different ligands which bind to the ATP-binding site differently modulate responses of allosteric site in interaction with a partner protein (e.g., enhance or inhibit interaction). Providing computer platform for rational design of allosteric modulators, the authors of this study opened an exciting avenue for future discovery of novel class of kinase inhibitors with less on-target side effects and more specific modulation of signaling pathways.

#### Case Study

#### **Rational design of clinical candidate Asciminib—allosteric Bcr-Abl1 inhibitor**

Asciminib belongs to a class of drugs designed to inhibit Bcr-Abl by binding to an allosteric pocket known as myristatebinding pocket. Rational development of Asciminib started with fragment-based screening using NMR assay (Schoepfer et al., 2018). Although determined NMR-based dissociation constants (Kd) for fragment hits were satisfactory, none of the fragments were active in biochemical and cellular assays. Subsequent X-ray studies revealed inability of fragment hits to induce assembled inactive state by bending of helix I, previously reported as conformational change important or autoinhibition of Abl by myristoilation (Nagar et al., 2003). Following this finding, the authors established another screening assay, the NMR-based conformational assay, which monitors the conformational state of C-terminal helix I (Jahnke et al., 2010). NMR-based conformational assay was used to investigate identified fragments and series of known allosteric modulators derivatives of GNF-2 (Adrian et al., 2006; **Figure 9**). Results of the study revealed that compounds which bind to myristoyl pocket and do not induce helix I bending were actually functional activators of Abl1 (by interfering with autoinhibition mechanism of Abl1). Critical bending of helix I was found to be induced by the presence of CF3O– group from GNF-2. Based on these findings, CADD techniques (molecular docking, similarity and pharmacophore searches) were used to design compound **X** in respect to X-ray structure with bent helix I conformation. Subsequent introduction of CF3O– group finally led to the first active allosteric inhibitor. Molecular modeling techniques were used in combination with X-ray crystallography in order to optimize potency and drug-like properties of the compound. Although only standard CADD techniques were reported in the discovery of Asciminib, recent application of molecular dynamics-based approaches demonstrated utility of such techniques in examination of mechanisms of resistance

5 and finally clinical candidate—Asciminib.

and effects of dual targeting of ATP-binding and allosteric site providing rationale for development of novel drugs (El Rashedy et al., 2018; Meng et al., 2018; Zhan et al., 2019).

#### **Machine learning methods to predict kinase-compound interactions**

Nowadays, we are seeing the widespread use of machine learning in many areas, including pharmaceutical industry, especially in drug design. Popular computational methods initially used in pharmaceutical research were quantitative structure activity relationships (QSAR) and quantitative structure property relationship (QSPR), which were adequate for small datasets. However, with the rapid growth of databases (thanks to methods such as high-throughput in vitro screening and Xray crystallography), it became inevitable to develop different in silico tools that can manage bigger data (Ekins, 2016). Today, many different machine learning methods such as support vector machines (SVM), k-Nearest Neighbors, Artificial Neural Networks (ANN), Deep Learning (DL), etc. are used in pharmaceutical research and they can be applied in various processes of drug design from virtual screening to de novo drug design (Buchwald et al., 2011; Drewry et al., 2017; Konze et al., 2019; Kuthuru et al., 2019; Lee et al., 2019; Zhavoronkov et al., 2019).

Many different machine learning models were created for the prediction of drug–target interactions (DTI), and many DTI methods have been applied to the protein kinases family (Kuthuru et al., 2019). Unlike LB and SB methods, DTI prediction uses the information from both protein and ligand and these methods can be similarity based or descriptor-based. One of the first similarity-based methods for identification of drug–target interactions was introduced by Yamanishi et al. in 2008. It used the known drug structure, protein sequence and drug–target interaction network to determine unknown ligand–target interactions. The main hypothesis is that two compounds that have high structure similarity might probably interact with similar target proteins, and vice versa two proteins with high sequence similarity might probably interact with similar drugs (Yamanishi et al., 2008). On the other hand, descriptor-based models use feature vectors from known drug structures and protein sequences, as inputs for machine learning methods, such SVM, AAN, DL, etc. In 2011, Buchwald et al. used SVM to prepare the model for prediction of protein kinases–ligand interactions. They used a set of binding data obtained from 113 different protein kinases and 20 inhibitors obtained through ATP site-dependent binding competition assays. They focused on vector features that describe the structure of molecules that are connected with certain chemical environment–protein active site sequence and created a SVM model with good predictivity (Buchwald et al., 2011).

Recently, the use of ANN, especially deep learning methods saw a significant increase in the process of drug design (Ekins, 2016; Merk et al., 2018; Putin et al., 2018; Konze et al., 2019). Deep generative models are utilizing neural networks to generate new objects (drugs) with desired properties (for example activity, Ki, IC50). These methods should be able to produce chemically correct structures without the need for including fragment libraries and/or rules for their combination (Merk et al., 2018). The ability to produce novel chemical structures with certain properties makes deep generative models suitable for the discovery of novel possible therapeutics (Zhavoronkov et al., 2019). In 2018, Merk et al. applied generative models to come up with novel bioactive, synthesizable drugs. They trained the model with more than 500,000 SMILES of bioactive compounds with their activity properties extracted from the ChEMBL (KD, Ki , IC/EC50 values <1µM). Additionally, the model was fine-tuned to enable the de novo generation of targetspecific ligands on retinoid X receptors (RXR) and/or peroxisome proliferator-activated receptors (PPAR). Finally, none of the generated compounds was identical to compounds from the training sets, and they were residing within the RXR/PPAR region of the fine-tuning set (Merk et al., 2018).

#### Case Studies

#### **Predictive proteochemometric models for kinases derived from 3D protein field-based descriptors**

Subramanian et al. described the development of proteochemometric models for 1,572 inhibitors and 95 kinases obtained from Kinase SARfari (https://chembl.gitbook. io/chembl-interface-documentation/legacy-resources#kinase-

sarfari) and CHEMBL database, using 3D structure of proteins and active and inactive ligands. Proteins were described with molecular interaction fields derived from Schrödinger's WaterMaps, while different 1D, 2D, and 3D descriptors were used to describe the ligands. Separate training sets were created for the ligands and targets. Different methods were used for preparation of the proteochemometric models: support vector machines (SVM) and random forests (RF). The ligand prediction model was trained on the ligand training set and was used for ligand prediction model and target training set for target predicting model. In the end, they validated all the models using internal and external validation. This approach allows creation of not only predictive proteochemometrics model for protein kinases, but also preparation of visually interpretable models. This allows interpretation of kinase–ligand interactions, which can be used, for example, for optimization of ligand in order to achieve optimal activity and/or selectivity. Having visually interpretable models is the advantage compared to classical DTI methods that use only 2D information (Subramanian et al., 2013, 2016).

#### **Deep learning model for identification of potent discoidin domain receptor 1 kinase inhibitors**

Recently, Zhavoronkov et al. created a deep generative model for de novo small-molecule design—GENTRL (GENerative Tensorial Reinforcement Learning). Besides the effectiveness of a compound against a given biological target, GENTRL also takes into account its dissimilarity from other molecules in the literature and patent space, as well as its synthetic feasibility. For the proof-of-concept GENTRL was used to design potential Discoidin domain receptor 1 (DDR1) kinase inhibitors. Data was collected from different data sets: ZINC data set, known DDR1 kinase inhibitors data set, common kinase inhibitors, molecules with activity on non-kinase targets, patent data, and used to train the model. The model was generated by combining reinforcement learning with a reward, variational inference, and tensor decompositions. Finally the randomly elected six compounds that have not been previously published or patented were designed, synthesized, and experimentally tested. The whole process lasted only 46 days, which suggests that the application of drug design methods such as this will reduce the time and cost of drug discovery process (Zhavoronkov et al., 2019).

#### CONCLUDING REMARKS

In silico approaches are viable, usually cheaper and faster alternative to experimental drug discovery techniques. This review summarizes the most important computational tools that have led to the discovery of kinase inhibitors, many of which are in clinical use today as promising anticancer drugs. Computational approaches, such as QSAR modeling, ligandbased and structure-based virtual screening, molecular docking, molecular dynamics, quantum mechanics, fragment-based drug design, and machine learning methods, provide unique insight in the conformational landscape of kinases, structural requirements for inhibitory activity, binding modes and atomistic mechanisms of allostery, which represent indispensable information for rational de novo design. One of the main advantages of computational approaches is the possibility of introduction of new groups on the known scaffolds and in silico prediction of activities and binding affinities. Known scaffolds of the approved KIs include pyrimidine (imatinib, dasatinib, nilotinib), quinazoline (erlotinib, gefitinib, afatinib, vandetanib), pyridine (sorafenib), pirrolopyridine (vemurafenib), pyrazolopyridine (ibrutinib) etc. In silico modification of these scaffolds resulted in the design of many kinase inhibitors with enhanced predicted activities and binding affinities which can serve as lead compounds for further synthesis and preclinical testing. New chemical scaffolds that possess kinase inhibitory activity (imidazopyridazine, imidazopyridine, isoquinoline, phenazinamine, etc.) have also been proposed by computational approach and represent a good starting point for discovery of new kinase inhibitors. Due to increases in computational power, algorithmic improvements and increased accuracy, in silico approaches are yet expected to radically shape the era of kinase inhibitor discovery. Of note is to emphasize that not all drug discovery projects could be initiated and guided only with computational studies. The computational chemist must be aware of the structural biology of the studied targets, their dynamical changes influenced upon fragment/ligand binding. Whenever possible, it is advised to start CADD studies with experimental data and continue in silico optimization with combined modeling approaches, as much as possible. This review highlights the recent advances in discovery of kinase inhibitors by in silico approaches and can be useful for future design and synthesis of new kinase inhibitors as anticancer drugs.

#### AUTHOR CONTRIBUTIONS

ZG wrote the introduction and LB methods. TD wrote the structure of protein kinases, machine learning methods, and SB methods. DR wrote FB methods. ND wrote modeling of

#### REFERENCES

allosteric kinase inhibitors. KN contributed to the conception of the manuscript, collection of the data, drafting and revising. All authors have contributed to the bibliographical research and interpretation of the work, to its critical revision, and approved the final version of the manuscript.

#### ACKNOWLEDGMENTS

The authors kindly acknowledge national project number 172033 supported by the Ministry of Education, Science and Technological development of the Republic of Serbia.

hybrids as multipotent cholinesterase/monoamine oxidase inhibitors for the potential treatment of Alzheimer's disease. Eur. J. Med. Chem. 75, 82–95. doi: 10.1016/j.ejmech.2013.12.028

inhibitor set enabling the accurate prediction of compound-kinase activities: a way toward selective promiscuity by design?J. Chem. Inf. Model. 56, 1654–1675. doi: 10.1021/acs.jcim.6b00122

inhibitors: a rational approach to anticancer drug design. J. Saudi Chem. Soc. 17, 361–379. doi: 10.1016/j.jscs.2011.04.017

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Gagic, Ruzic, Djokovic, Djikic and Nikolic. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Vienna LiverTox Workspace—A Set of Machine Learning Models for Prediction of Interactions Profiles of Small Molecules With Transporters Relevant for Regulatory Agencies

Floriane Montanari † , Bernhard Knasmüller, Stefan Kohlbacher, Christoph Hillisch, Christine Baierová, Melanie Grandits\* and Gerhard F. Ecker

Pharmacoinformatics Research Group, Department of Pharmaceutical Chemistry, University of Vienna, Vienna, Austria

#### Edited by:

Jose L. Medina-Franco, National Autonomous University of Mexico, Mexico

#### Reviewed by:

Vinicius M. Alves, University of North Carolina at Chapel Hill, United States Rodolpho C. Braga, InsilicAll, Brazil

#### \*Correspondence:

Melanie Grandits<EMAIL_ADDRESS>
#### †Present address:

Floriane Montanari, Department of Digital Technologies, Bayer AG, Berlin, Germany

#### Specialty section:

This article was submitted to Medicinal and Pharmaceutical Chemistry, a section of the journal Frontiers in Chemistry

Received: 09 October 2019 Accepted: 13 December 2019 Published: 10 January 2020

#### Citation:

Montanari F, Knasmüller B, Kohlbacher S, Hillisch C, Baierová C, Grandits M and Ecker GF (2020) Vienna LiverTox Workspace—A Set of Machine Learning Models for Prediction of Interactions Profiles of Small Molecules With Transporters Relevant for Regulatory Agencies. Front. Chem. 7:899. doi: 10.3389/fchem.2019.00899 Transporters expressed in the liver play a major role in drug pharmacokinetics and are a key component of the physiological bile flow. Inhibition of these transporters may lead to drug-drug interactions or even drug-induced liver injury. Therefore, predicting the interaction profile of small molecules with transporters expressed in the liver may help medicinal chemists and toxicologists to prioritize compounds in an early phase of the drug development process. Based on a comprehensive analysis of the data available in the public domain, we developed a set of classification models which allow to predict—for a small molecule—the inhibition of and transport by a set of liver transporters considered to be relevant by FDA, EMA, and the Japanese regulatory agency. The models were validated by cross-validation and external test sets and comprise cross validated balanced accuracies in the range of 0.64–0.88. Finally, models were implemented as an easy to use web-service which is freely available at https:// livertox.univie.ac.at.

Keywords: Vienna LiverTox Workspace, web service, machine learning, ABC-transporter, OATP-transporter, toxicity, classification models

### INTRODUCTION

Membrane transporters expressed in the liver play different, but interconnected roles: on the one hand, basolateral transporters pick up xenobiotics, and endogenous molecules from the portal vein to the liver, or excrete their substrates into the blood. Apical transporters, on the other hand, take care of the flux toward the bile duct network (**Figure 1**). Three main types of substrates are of interest with respect to liver toxicity: drugs, which enter the hepatocytes at the first hepatic pass or at the elimination stage; bilirubin, a product of the degradation of the heme; and bile salts, which circulate between the gastro-intestinal tract, and the liver.

Additionally, to the enzyme family of cytochromes, also the transporters expressed in the liver are crucial for a fully functional organ. Some of them are e.g., involved in the bilirubin cycle: OATP1B1 and OATP1B3 uptake bilirubin into the hepatocytes (Briz et al., 2003), where glucuronidation takes place. MRP2 then excretes the bilirubin conjugate to the bile (Kamisako et al., 1999). At the basolateral membrane, MRP3 might also excrete it back to the sinusoidal blood (Keppler, 2014). As a result, inhibition of the uptake OATP transporters or of MRP2 may

lead to the accumulation of bilirubin (conjugated or not) in the blood, which is referred to as hyperbilirubinemia. Conjugated hyperbilirubinemia is a marker of hepatobiliary injury (Dufour et al., 2000; Ozer et al., 2008; Padda et al., 2011), and predicting it may allow to flag compounds that could cause liver injury.

Bile acids are synthesized in the liver by catabolism of cholesterol and then excreted to the bile by the active bile salt export pump (BSEP) and by the multidrug resistance-associated protein 2 (MRP2) (Meier and Stieger, 2002). Bile salts have a pronounced detergent effect, which explains their toxicity when they accumulate in the liver (Attili et al., 1986). For their transport in the bile duct, bile salts form mixed micelles with phospholipids of the outer leaflet of the membrane. The multidrug resistance protein 3 (MDR3) allows lipid flopping at the apical membrane of the hepatocyte, and its function is necessary to avoid bile duct toxicity (Nicolaou et al., 2012). After reaching the intestine via the bile flow, bile acids are reabsorbed into the portal vein, and taken up again into the hepatocytes by the sodium taurocholate co-transporting polypeptide (NTCP) (Stieger, 2011). Impairment of bile flow leading to a toxic accumulation of bile salts in the hepatocytes might lead to druginduced cholestasis, which is one of the main causes of druginduced liver injury (DILI) (Padda et al., 2011).

Apart from playing a role in proper bile flow and bilirubin elimination, liver transporters also transport drugs that will then be metabolized and excreted. At this stage, drugs can inhibit different transporters and cause drug-drug interactions (König et al., 2013) (in case of co-administered inhibitor and substrate) or liver injury (by disrupting the bile flow for example). This is why predicting the inhibition and the substrate profile for liver transporters might be useful in identifying potentially problematic compounds. In addition, the Food and Drugs administration recommends experimental testing of the interactions between drugs and transporters (especially P-gp, BCRP, OATP1B1, and OATP1B3) to identify potential drugdrug interactions (U.S. Department of Health Human Services Food Drug Administration Center for Drug Evaluation Research, 2012). Thus, it definitely would be of value to have a suite of computational models available which allow the fast and easy assessment of compounds for their interaction profile with transporters expressed in the liver.

Here, we present the Vienna LiverTox Workspace, a web server for the prediction of interactions with liver transporters as well as selected liver toxicity endpoints. To the best of our knowledge, it is the first time that such an ensemble of predictive models for hepatotoxicity and liver transport is made available to the public. The predictions are made by individual machine learning models built on publicly available data for each target of interest.

### METHODS

### Data Curation

The datasets for the training and testing of the models were collected from different sources (online tools as well as publications). The data were cleaned by using an in-house system combining the molecular operating environment (MOE 2014.09) (Molecular Operating Environment, 2014) wash option and the Atkinson Standardiser (https://github.com/flatkinson/ standardiser). This approach was used for some of the datasets and for others the cleaning procedure of already published papers were used (Pinto et al., 2012; Kotsampasakou et al., 2015). In general, duplicates were removed from the dataset, including pairs of stereoisomers. Further, if these compounds share the same class label, one of the compounds was kept. A detailed list of references as well as the number of compounds revealed after the preprocessing of the data is given in **Table S1**. In the sections Transporter Models and Hepatotoxicity Models, the data curation as well as the model generation for the specific endpoints is given. The datasets are available in the Supplementary Material (**Data Sheet 2**).

#### Transporter Models

The web service allows for the prediction of interactions between a small molecule and eight different liver transporters (**Figure 1**, transporters marked in blue). The lack of publicly available data for the other transporters explains the absence of respective models in the Vienna LiverTox Workspace.

All the models predicting whether a compound will be a substrate of a transporter (BCRP, P-gp, BSEP, MRP2, and MRP3) were built upon a dataset correlating expression levels of 47 ABC-transporters with drug toxicity, which can then be used to infer transported vs. non-transported compounds (Szakács et al., 2004). For the transport inhibition models (BCRP, P-gp, BSEP, MRP3, MRP4, OATP1B1, and OATP1B3) the datasets were collected from literature and, if necessary, manually aggregated. In both cases, the models predict a binary outcome: the query compound is a substrate or not, or an inhibitor or not.

For chemistry encoding of the compounds, we used circular fingerprints or 2D molecular descriptors as implemented in

**Abbreviations:** ADME-Tox, absorption distribution metabolism excretion toxicity; BCRP, breast cancer resistance protein; BSEP, bile salt export pump; DILI, drug-induced liver injury; ECFP, extended connectivity fingerprint, MDR1, multidrug resistance protein 1; MDR3, multidrug resistance protein 3; MRP2, multidrug resistance-associated protein 2; MRP3, multidrug resistance-associated protein 3; MRP4, multidrug resistance-associated protein 4; NTCP, sodium taurocholate co-transporting polypeptide; OATP1B1, organic anion transporting polypeptide 1B1; OATP1B3, organic anion transporting polypeptide 1B3; P-gp, P-glycoprotein; SMILES, simplified molecular input line entry system.

\*see https://livertox.univie.ac.at/ for a detailed list.

RDKit version 2015.03.1 (https://www.rdkit.org/). Different machine learning algorithms were applied and the one giving the best cross-validation results was kept as final model. Especially for the transport prediction, a heavy class imbalancy (most of the drugs in the training set were non-substrates) was noted, which was handled by MetaCost (Domingos, 1999). The exact methodology and cross-validation performance for each individual transporter model is described in the documentation available at https://livertox.univie.ac.at, and an overview is given in **Table 1**. In some cases, external test sets were collected from Metrabase (Mak et al., 2015) or from recent publications (timesplit evaluation).

#### Hepatotoxicity Models

Three models in the web service can be used to assess human liver damage potentially caused by a test compound: hyperbilirubinemia, cholestasis, and drug-induced liver injury (DILI).

For hyperbilirubinemia, 835 compounds were taken from Kotsampasakou et al. (2017a) and the modeling methodology was kept as in Kotsampasakou et al. (2017b): ECFP-like fingerprints were computed with RDKit (https://www.rdkit. org/), then a combination of feature selection, MetaCost, and support vector machines with RBF kernel was used for learning. The cholestasis model uses data from Kotsampasakou and Ecker (2017) and a combination of MetaCost and a tree algorithm to predict whether a compound is likely to cause cholestasis or not. Finally, the DILI model is based on a 966-compound dataset carefully compiled from literature (Kotsampasakou et al., 2017c). RDKit molecular descriptors and a random forest of 500 trees are used for modeling.

#### Web Service Implementation

The Vienna LiverTox Workspace has been implemented as Python/php based web service. It consists of two parts, namely the backend and the frontend. The backend consists of a docker image which runs the machine learning models on an input SD-File. It consists of a Python Flask server (https://palletsprojects. com/p/flask/) which processes the requests from the frontend. Each request consists of one or more input molecules and a list of models to run the predictions on. The frontend, also a docker image, is based on the CakePHP framework (https://cakephp. org/) and is responsible for the user interface (UI), which sends the request to the backend and displays the results. The web service provides, after a login, the possibility to upload a SD-File of 10 compounds. The service can also be used without logging in but then it is only possible to draw and predict a single molecule. JSME (Bienfait and Ertl, 2013) is used as drawing tool. The web service runs on an Ubuntu Linux based server with two twelve-core Intel Xeon 64bit processors and 128 GB RAM, and is hosted at the University of Vienna by the Pharmacoinformatics Research Group.

The models use the RDKit library (https://www.rdkit.org/) (version 2015.03.1) for computing the descriptors and handling the chemistry aspects, while the Weka (Hall et al., 2009) (version 3.7.11) and scikit-learn (Pedregosa et al., 2011) (version 0.14.1) libraries are used to train and run the predictive models. The models also include a compound cleaning step, implemented with the Atkinson Standardiser (https://github.com/flatkinson/ standardiser) (**Figure 2**).

The steps performed during model building and a test compound prediction are shown in **Figure 2**. In both cases, the compound is standardized and the molecular descriptors are calculated. In the case of the model generation, this allows the training of the model and its development (left workflow). For the prediction of a test compound, the descriptors are passed to the available model to predict its class affiliation (right workflow).

In general, the output of the model gives, in addition to the class prediction, the actual score. This score is a numerical value between 0 and 1, and roughly corresponds to a probability of being active (inhibitor, substrate, or toxic compound). Therefore, a value close to 1 indicates substrate/inhibitor/toxic properties, a value close to 0 annotates for non-substrates/noninhibitors/non-toxic.

#### Applicability Domain

The Applicability Domain (AD) is used to validate the reliability of a given prediction model. It defines whether

TABLE 2 | Performance metrics of the transporter models.

a dataset of interest is in or out of domain, meaning if it falls within the chemical space of the model or not. If it is out of Domain, the prediction cannot be regarded as reliable.

In our study, an Applicability Domain model, using the approach of Sahigara et al. (2013), was created for each transporter with the respective training set. RDKit descriptors were used as molecule representation. This approach combines the classical, widely used k-Nearest Neighbor (k-NN) method with a probability density function estimation. It uses three stages to determine the reliable space of a prediction model. First, a set of thresholds is defined depending on the diverse densities of the training set by considering the 15 k-nearest neighbors using Euclidean distances. This allows the AD to consider a dense and sparse training region (The threshold defines if a test sample can be reliable predicted). In a next step a decision rule is derived to filter out outlier molecules. Finally, the reliability of the AD is tested by looking at the model statistics and prediction errors. This feature is not yet implemented in the Web service, but will be available soon.

### RESULTS

#### Model Results

The performance of the models was estimated by calculating statistical performance metrics using a 10-fold cross-validation. The results are provided in **Table 2**. The overall accuracy, corresponding to the rate of correct predictions, ranges from 0.59 to 0.87. Also, the sensitivity of the models was calculated (0.57–0.85). This parameter gives the number of actual positives that are correctly identified and is expressed by the number of true positives divided by the number of positive predictions. Further, the number of actual negatives was determined by the

FIGURE 3 | Overview of the web service interface. On the left side, test compound is drawn and desired models are selected. On the right side, results table with the predictions and model scores.

number of true negatives divided by the number of negative predictions. The so-called specificity ranges from 0.56 to 0.90. To estimate a metric for the quality of the models, the Matthews correlation coefficient (MCC) and the Area under the Receiver Operating Characteristics curve (ROC AUC) were determined. The MCC is a number between −1 and 1 where 0 indicates a prediction equal to a random prediction and 1 indicates a perfect prediction, whereas −1 is a complete miss. The scores for our models are in the range of 0.20–0.76. The ROC AUC measures the ability of the model to distinguish between negatives and positives, while a higher value indicates a better performance. In the models provided on the Vienna LiverTox Workspace the ROC AUC ranges from 0.64 to 0.94. Furthermore, if data was available, the models were also validated with one or more external test sets. For more details see the documentation on the website.

### Use Case: Prediction of Liver Interaction for a Propafenone Analog

In this section, we briefly detail how predictions can be generated for a given compound using the LiverTox web service. In first instance, the compound is drawn or its SMILES string is pasted in the Molecule Editor (Bienfait and Ertl, 2013). Then the models can be selected on the left panel, either one by one or all at the same time (**Figure 3**, left side).

By clicking on the SUBMIT button, the data is sent to the backend server where the predictions are running. Upon completion of the calculations, a table listing the different models and corresponding outputs will appear (**Figure 3**, right side). The second column in the results table corresponds to the binary classification, while the third column "Score" gives the actual output of the model, which corresponds to a probability of being active. For example, for BCRP inhibition and transport TABLE 3 | Comparison of existing free online tools to predict ADME-Tox properties of compounds.

and for DILI, the output score is close to 0.5 (which is the threshold used to separate predicted actives from predicted inactives), which indicates an uncertainty of these three models for the particular query compound. Propafenone derivatives are frequently reported as inhibitors of P-gp, and indeed the P-gp inhibition model predicts this particular one to be an inhibitor with a high score (0.93).

#### DISCUSSION

Many systems already exist to predict in silico activities or properties of small molecules. **Table 3** compares freely available ones with our own web service in terms of model offer, submission and run time. For example, ProTox-II predicts oral drug toxicity in rodents (lethal dose LD50 and a category of toxicity between 1 and 6) using similarity to compounds with known LD50 and recognition of toxic fragments (Drwal et al., 2014). BioZyne proposes exclusively one model for P-gp transport prediction based on the same dataset as ours (Szakács et al., 2004; Levatic et al., 2013 ´ ). It uses a Support Vector Machine classifier for the prediction of P-gp substrates. The Danish (Q)SAR Database contains pre-calculated properties combined from more than 200 models from both commercial and free tools (http://qsar.food.dtu.dk/). Predictions for environmental toxicity, blood-brain barrier permeation, cytochrome interactions, or human genotoxicity are available. Unfortunately, new predictions for compounds that are not part of the database cannot be made. PkCSM is another web service for predicting pharmacokinetics properties of compounds (Pires et al., 2015). Models such as P-gp inhibition and transport, bloodbrain barrier permeation, interaction with cytochromes, renal clearance, or even liver toxicity are available.

In general, our models for the inhibitors show a better performance especially when looking at the correct prediction of the positives. The prediction of true negatives is for the inhibitor and transporter models quite similar which can be explained by the availability of more negatives if the training set is unbalanced. This is especially the case for the substrate models. The quality of the prediction (MCC) is higher for the inhibition models of P-gp, BSEP, BCRP, and MRP3 since the available dataset is more balanced. In comparison, the three toxicity models show a poorer performance due to the complexity of these endpoints and especially for hyperbilirubinemia and cholestasis which shows also a lack of positives.

The Transporters selected for this web service were chosen based on their importance for regulatory agencies such as FDA, EMA and the Japanese regulatory agency. They recommend or in some cases request these proteins to be routinely tested in inhibition—and substrate studies of new drugs.

#### CONCLUSION

We have presented the Vienna LiverTox Workspace, a web service dedicated to the prediction of liver toxicity and interactions between small molecules and liver transporters. It is easy to use, fast, web browser agnostic, and well-documented. Thanks to its modular system, it will be easy to integrate new models in the future, as well as re-implement existing models in case new training data becomes available. We hope that our models will help researchers to flag potentially dangerous compounds and shed light on the relationships between liver transporters and toxicity.

#### DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/**Supplementary Material**.

#### AUTHOR CONTRIBUTIONS

FM, SK, CH, and MG developed the models. BK, CB, and MG implemented the web service. MG and GE supervised the study. FM and MG wrote the majority of the manuscript. All authors contributed to refining the manuscript.

#### FUNDING

We gratefully acknowledge financial support provided by the Austrian Science Fund, grant #F03502 (SFB35), and by the Innovative Medicines Initiative Joint Undertaking under grant agreement n◦ 115002 (eTOX).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2019.00899/full#supplementary-material

### REFERENCES

profiling ABC transporter genes in cancer cells. Cancer Cell 6, 129–37. doi: 10.1016/j.ccr.2004.06.026

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Montanari, Knasmüller, Kohlbacher, Hillisch, Baierová, Grandits and Ecker. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# ChEMBL-Likeness Score and Database GDBChEMBL

#### Sven Bühlmann and Jean-Louis Reymond\*

Department of Chemistry and Biochemistry, University of Bern, Bern, Switzerland

The generated database GDB17 enumerates 166.4 billion molecules up to 17 atoms of C, N, O, S and halogens following simple rules of chemical stability and synthetic feasibility. However, most molecules in GDB17 are too complex to be considered for chemical synthesis. To address this limitation, we report GDBChEMBL as a subset of GDB17 featuring 10 million molecules selected according to a ChEMBL-likeness score (CLscore) calculated from the frequency of occurrence of circular substructures in ChEMBL, followed by uniform sampling across molecular size, stereocenters and heteroatoms. Compared to the previously reported subsets FDB17 and GDBMedChem selected from GDB17 by fragment-likeness, respectively, medicinal chemistry criteria, our new subset features molecules with higher synthetic accessibility and possibly bioactivity yet retains a broad and continuous coverage of chemical space typical of the entire GDB17. GDBChEMBL is accessible at http://gdb.unibe.ch for download and for browsing using an interactive chemical space map at http://faerun.gdb.tools.

#### Edited by:

Jose L. Medina-Franco, National Autonomous University of Mexico, Mexico

#### Reviewed by:

Cristian G. Bologa, University of New Mexico, United States Rodrigo Ochoa, University of Antioquia, Colombia

#### \*Correspondence:

Jean-Louis Reymond<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Medicinal and Pharmaceutical Chemistry, a section of the journal Frontiers in Chemistry

Received: 25 October 2019 Accepted: 15 January 2020 Published: 04 February 2020

#### Citation:

Bühlmann S and Reymond J-L (2020) ChEMBL-Likeness Score and Database GDBChEMBL. Front. Chem. 8:46. doi: 10.3389/fchem.2020.00046 Keywords: chemical space exploration, molecular database, enumeration algorithm, chemical space mapping, virtual screening

#### INTRODUCTION

Innovation at the level of chemical structures is an essential part of drug discovery. Novelty often results from chemical intuition however this approach is increasingly difficult as the number of known molecules increases. Novelty is similarly limited in virtual combinatorial libraries (Leach and Hann, 2000; Hu et al., 2011; van Hilten et al., 2019) and generative models trained with known molecules (Chen et al., 2018; Elton et al., 2019) because these systems mostly shuffle known patterns, which produces many technically new but often not fundamentally innovative molecules. To circumvent this limitation, we have initiated the exhaustive enumeration of all possible organic molecules following simple rules of chemical stability and synthetic feasibility, and reported large databases enumerating molecules up to 11 (Fink et al., 2005; Fink and Reymond, 2007), 13 (Blum and Reymond, 2009), and 17 atoms (Ruddigkeit et al., 2012, 2013), as well as of possible ring systems up to 30 atoms (Visini et al., 2017a). Analyzing the resulting generated databases (GDBs) shows that there are many orders of magnitude more possible molecules spanning a much broader structural diversity than already known ones (Reymond, 2015; Awale et al., 2017b).

One of the defining features of the GDB databases is the exponential increase in the number of possible molecules as function of increasing molecular size and complexity elements, such as stereocenters and heteroatoms, implying that most possible molecules are in fact far too complex to be considered as realistic synthetic targets. To address this problem we have designed subsets of our largest database GDB17 by limiting complexity elements using simplification criteria, such as fragment-likeness (Congreve et al., 2003), producing the fragment database FDB17, and

**290**

medicinal chemistry rules for functional groups and complexity (Mignani et al., 2018), producing the medicinal chemistry aware database GDBMedChem (Visini et al., 2017b; Awale et al., 2019). These approaches however also constrain the diversity of GDB molecules, which partly defeats the purpose of exploring chemical space broadly.

Herein we report an alternative approach to create subsets of GDB17 based on the frequency of occurrence of substructures from known molecules independent of the overall molecular structure (**Figure 1A**). We define a "ChEMBL-likeness" score (CLscore) by considering which substructures in a molecule also occur in molecules from the public database ChEMBL (Gaulton et al., 2017), using a subset of molecules with reported high confidence datapoint of activity on single protein targets, a type of ChEMBL subset which we have used previously for target prediction (Awale and Reymond, 2019; Poirier et al., 2019). We then filter the entire GDB17 with a cut-off value for CLscore, followed by uniform sampling of the resulting subset across molecular size, stereocenters and heteroatoms as done previously with FDB17 and GDBMedChem, to obtain a ChEMBL-like subset of 10 million molecules forming the database GDBChEMBL. This database covers chemical space as broadly as but more continuously than FDB17 and GDBMedChem yet features a much higher synthetic accessibility as judged by a calculated synthetic accessibility score (Ertl and Schuffenhauer, 2009), might contain molecules with a higher probability of bioactivity, and in any case provides a very different starting point to serve as a source of inspiration for molecular design.

#### RESULTS AND DISCUSSION

#### ChEMBL-Likeness Score

Our definition of CLscore is related to the synthetic accessibility score (SAscore) (Ertl and Schuffenhauer, 2009) and natural product likeness score (NPscore) (Jayaseelan et al., 2012) of a molecule, which are calculated from the occurrence frequencies of its substructures in PubChem and fragments from natural products, respectively, combined with additional functional group rules. Here we focus on 457,139 compounds recorded in ChEMBL24 as being active on single protein targets (IC50 or EC50 ≤ 10µM) with high confidence datapoints (Awale and Reymond, 2019; Poirier et al., 2019). To design our CLscore we consider circular substructures, called molecular shingles, because they form the basis for molecular fingerprints ECFP4 and MHFP6 which perform best in benchmarking studies (Riniker and Landrum, 2013; Probst and Reymond, 2018).

The frequency of occurrence of the 636,979 molecular shingles up to a diameter of six bonds found in our ChEMBL subset follows a power law distribution (**Figure 1C**). To compute the CLscore of a molecule, we assign to each of its shingles (S) a shingle value calculated from the logarithm of its frequency of occurrence fS in our ChEMBL subset, considering only shingles occurring at least 100 times in this subset (141,261 shingles, 22.2% of the total). We then sum all shingle values and divide the sum by the total number of shingles in the molecule (Equation 1).

$$\text{CLscore} = \frac{\sum\_{i=1}^{m} \log\_{10} \left( f\_{\mathcal{S}} \right)\_{i}}{\text{N}}\\\begin{array}{l} \text{CLscore} := \text{ChEBML} - \text{likeness score} \\ \text{S} := \text{shningleinmolecular structure} \\ f\_{\mathcal{S}} := \text{abundanceofmolecularishingleinChEMBL} \\ \text{N} := \text{totalnumberofmightsholder structure} \\ \text{m} := \text{numberofmingslandmeedwith ChEMBL} \end{array}$$

The histogram of CLscore for the 457,139 ChEMBL reference molecules is approximately Gaussian with a peak at CLscore = 3.9 (**Figure 1B**). DrugBank (Law et al., 2014) and particularly ZINC (Sterling and Irwin, 2015) peak at a similar CLscore, showing that these three databases consist of molecules built from the same type of substructures. By contrast GDB17 and its subsets FDB17 and GDBMedChem have a much lower CLscore distribution peaking at CLscore = 2.7, reflecting the fact that GDB molecules are very different from ChEMBL molecules. CLscore values correlate with SAscore values, reflecting the similar principles underlying both scores, and suggesting that molecules with high CLscore should also be synthetically accessible (**Figure 1D**).

#### GDBChEMBL Database

Calculating CLscores on the entire GDB17 (166.4 billion SMILES) and keeping molecules with CLscore ≥ 3.3, a cut-off value which retains 78.3% of our ChEMBL subset, eliminates 84.3% of GDB17. The remaining 26.2 billion molecules are then binned in triplet value bins considering heavy atom count (HAC 1-17), stereocenter count (0–4, ≥ 5) and heteroatom count (0– 8, ≥ 8). There are 538 different triplet value bins, which are occupied by 1 to 1.6 × 109 molecules. Uniform sampling finally yields a final set of 10 million molecules evenly distributed across molecular size, stereochemical complexity and polarity, forming the database GDBChEMBL (**Figure 1E**).

As a consequence of uniform sampling, the heavy atom count (HAC) profile of GDBChEMBL resembles that of FDB17 and GDBMedChem and is relatively flat compared to the very steep peak at HAC = 17 in the parent database GDB17 (**Figure 2A**). Uniform sampling also explains the rotatable bond count (RBC) profile in GDB subsets compared to GDB17 (**Figure 2B**), as well as the fact that the profiles of the three GDB subsets across these parameters are generally more similar to the profile of molecules up to 17 atoms in ChEMBL (ChEMBL17) and to natural products (UNPD17) (Banerjee et al., 2015) than to the profile of GDB17.

GDBChEMBL displays a very broad distribution in terms of hydrogen bond donor atoms (HBD, **Figure 2C**), hydrogen bond acceptor atoms (HBA, **Figure 2D**) and nitrogen plus oxygen atom count (N+O, **Figure 2E**) due to the absence of heteroatom capping criteria in selecting GDBChEMBL compared to FDB17 and GDBMedChem, for which fragment-likeness criteria, respectively, caps on the number of functional groups were applied. Similar differences are visible in topological polar surface area (TPSA, **Figure 2F**) and calculated octanol/water partition coefficient (alogP, **Figure 2G**). The broader distribution of polarity parameters in GDGChEMBL compared to GDB17 results from uniform sampling since the procedure gives relatively more importance to molecules with extreme size and polarity values.

Synthetic accessibility is better (lower SAscore) in GDBChEMBL than for GDB17, FDB17, or GDBMedChem,

reflecting the correlation between CLscore and SAscore noted above (**Figure 2H**). Similar to GDB17 and its other subsets, GDBChEMBL displays a much higher fraction of sp3 atoms than ChEMBL (fsp3 , **Figure 2I**). As a consequence GDB molecules are closer to natural products, which is reflected in the NPscore profile (**Figure 2J**). Despite of these differences and similarities in SAscore and NPscore, it must be noted that GDB17 and its subsets stand out by the fact that they contain fewer aromatic and more heterocyclic molecules than ChEMBL and natural products (**Figure 2K**).

#### Visualization and Similarity Searching

To gain an overview of GDBChEMBL we computed Molecular Quantum Number (MQN) fingerprint values (Nguyen et al., 2009), performed a principal component analysis (Rosén et al., 2009), and visualized the resulting 3D-map in the interactive web-based application faerun (Probst et al., 2018). In this 3Dmap accessible at http://faerun.gdb.tools, each point represents one or more molecules present at the corresponding position and can be color-coded according to a molecular property selected from the faerun menu.

Comparing MQN maps of GDBChEMBL, FDB17 and GDBMedChem shows that each of the three GDB17 subset cover a similar range of properties, however coverage by GDBChEMBL is more continuous, as is well visible in the vertical stripe at right containing all acyclic molecules (**Figures 3A–C**). Note that CLscore values are not correlated with MQN properties, which is not surprising considering that ChEMBL substructure span a broad range of properties (**Figure 3D**). Color-coding by the calculated logP value (alogP, **Figure 3E**) and by rotatable bond count (RBC, **Figure 3F**) illustrate the distribution of molecules in the MQN map.

The fact that molecules in GDBChEMBL are substantially different from those in the other subsets FDB17 and GDBMedChem can be shown by retrieving 1,000 MQN-nearest neighbors of nicotine from each database, and representing each dataset in a similarity map (Medina-Franco et al., 2007; Raghavendra and Maggiora, 2007; Awale and Reymond, 2015) using the molecular shape and pharmacophore fingerprint Xfp (Awale and Reymond, 2014), computed with the webbased application WebMolCS (Awale et al., 2017a). This visualization shows that each database provides different types of nicotine analogs (**Figure 3G**) with a good number of high

similarity analogs (**Figure 3H**). To facilitate similarity searches in GDBChEMBL, we have implemented a similarity search portal by which nearest neighbor searches of any molecule can be performed in GDBChEMBL using MQN, ECFP4, or a combined MQN-MHFP6 similarity, as described previously for GDBMedChem (Awale et al., 2019).

#### CONCLUSION

The data above demonstrate a substructure-based approach to select molecules from the generated database GDB17. As selection criterion we defined a ChEMBL-likeness score (CLscore) from the frequency occurrence of circular substructures, called molecular shingles, in a subset of the database ChEMBL consisting of compounds active on single protein targets with high confidence datapoints. This selection reduced GDB17 by 84.3%, leaving 26.2 billion molecules, which we sampled uniformly across molecular size, stereochemistry and heteroatoms to form GDBChEMBL comprising 10 million molecules.

Property profiles, chemical space maps and similarity searches show that GDBChEMBL is very different from our earlier GDB subsets FDB17 and GDBMedChem and spans chemical space more continuously. At the same time, the correlation between CLscore and the synthetic accessibility score (SAscore) implies that GDBChEMBL molecules will be on average easier

FIGURE 3 | Chemical space maps of GDBChEMBL, FDB17, and GDBMedChem. (A) PCA 3D-map of GDBChEMBL in MQN-space, color coded by heavy atom count; (B) same as a for FDB17; (C) same as a for GDBMedChem; (D) GDBChEMBL color-coded by CLscore value; (E) GDBChEMBL color-coded by calculated octanol/water partition coefficient alogP; (F) GDBChEMBL color-coded by rotatable bond count; (G) similarity map of MQN-nearest neighbors of nicotine from GDBChEMBL (red), FDB17 (cyan), and GDBMedChem (blue). Points in green and yellow indicate molecules shared by two databases. (H) Same as g color-coded by Xfp-similarity to nicotine. MQN maps a to f are accessible at http://faerun.gdb.tools. The similarity map of nicotine analogs g and h is accessible at: http://gdbtools. unibe.ch:8080/webMolCS/.

to synthesize than molecules from FDB17 and GDBMedChem, which have significantly lower CLscore and higher SAscores. We anticipate that the requirements for GDBChEMBL molecules to share a minimum number of substructures with molecules of known bioactivities from ChEMBL will also facilitate target prediction and the selection of interesting GDB molecules for synthesis and testing.

## METHODS

#### Preparative Steps ChEMBL Shingle Extraction

The ChEMBL (v 24.1) database was downloaded from https:// www.ebi.ac.uk/chembl/. Data points for extraction of molecular shingles were selected by applying the same restrictions that were used for extraction of training data for our Polypharmacology Browser PPB2 (Visini et al., 2017b). Structures were normalized to their major protonation state at pH 7.4 using ChemAxon cxcalc (v. 18.23.0). Molecular shingles for radii 1–3 were created using RDkit (2019.03.4) and converted to rooted, canonical, aromatic SMILES strings without retaining stereochemistry information. In association with abundancy in the ChEMBL, the SMILES substructures were stored as pickled python dictionary. Molecular substructures that were found <100 times were not stored.

#### CLscore Calculation

Scoring of GDB17 molecular structures was achieved by decomposition to molecular shingles in the exact same way as described for ChEMBL reference shingle extraction. For a specific query structure, all shingles are uniquely counted, then looked up in the ChEMBL reference database and upon match, logarithmic abundancy is summed up. The final CLscore is given by the ratio of total logarithmic abundancies of matched unique shingles to total unique shingles in the query structure. All respective scripts are accessible at: https://github.com/reymondgroup/GDBChEMBL.

#### GDBChEMBL Generation

All 166.4 billion molecular structures of GDB17 were decomposed to unique substructures in the same way as described for ChEMBL reference molecules. Only structures with CLscore ≥3.3 were stored. The final GDBChEMBL was obtained by distribution of all filtered 26.2 billion structures to 538 property triplet bins (heavy atom, heteroatom and stereocenter count). Property information was gathered using RDKit. Bins with 5+ hetero atoms and/or 8+ chiral atoms were merged. The actual even sampling was performed by sorting all property bins by size and defining target structure count as 10 million. Iteratively, remaining target count was divided by count of remaining bins, keeping all bins of size smaller than the current number to sample randomly. For each step, number of previously selected structures was subtracted from target count until random sample per remaining bins was lower than bin size. At this point, sample size was kept constant for all further bins.

#### REFERENCES

#### Visualizing GDBChEMBL in Faerun

Property color coded 3D maps for GDBChEMBL, FDB17, and GDBMedChem were generated using FUn (doc.gdb.tools/fun), an in-house developed framework for interactive visualization of chemical spaces on the web. Datasets were given as plain text, consisting of the four columns (space-separated): SMILES-string, numeric ID, 42 MQN descriptors (semicolon-separated) and further molecular properties used for map coloring (semicolonseparated). Next, the preprocessing toolchain was used to project the 42-dimensional MQN-space to 3D by applying Principal Component Analysis (PCA) and to generate all further files needed for visualization. Finally, the Underdark server was run using docker with Faerun visualization containers mapped.

#### Similarity Searching in GDBChEMBL

For better accessibility, GDBChEMBL is provided as a webbased interactive similarity search tool. The implementation uses HTML, Bootstrap, JavaScript, and the python Flask framework. Search times were reduced using Annoy trees (Approximate Nearest Neighbors Oh Yeah, https://github.com/spotify/annoy) which were created for the 42-dimensional MQN property space, as well as for 256-bit ECfp4. A third search option, MQN-MHFP6, initially searches using the MQN Annoy tree followed by resorting after Jaccard distance to query molecule in the MHFP6 fingerprint space (https://github.com/reymond-group/ mhfp). The search tool is available at: gdb.unibe.ch/tools.

#### DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/supplementary material.

#### AUTHOR CONTRIBUTIONS

SB designed and realized the study and wrote the paper. J-LR designed and supervised the study and wrote the paper.

#### FUNDING

This work was supported financially by the Swiss National Science Foundation.

as the target of a nanomolar angiogenesis inhibitor from a phenotypic screen using the polypharmacology browser PPB2. ChemMedChem 14, 224–236. doi: 10.1002/cmdc.201800554

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Bühlmann and Reymond. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Pilot Study of All-Computational Drug Design Protocol–From Structure Prediction to Interaction Analysis

#### Yifei Wu, Lei Lou and Zhong-Ru Xie\*

Computational Drug Discovery Laboratory, School of Electrical and Computer Engineering, College of Engineering, University of Georgia, Athens, GA, United States

Speeding up the drug discovery process is of great significance. To achieve that, high-efficiency methods should be exploited. The conventional wet-bench methods hardly meet the high-speed demand due to time-consuming experiments. Conversely, in silico approaches are much more efficient for drug discovery and design. However, in silico approaches usually serve as a supportive role in research processes. To fully exert the strength of computational methods, we propose a protocol which integrates various in silico approaches, from de novo protein structure prediction to ligand-protein interaction simulation. As a proof of concept, human SK2/calmodulin complex was used as a target for validation. First, we obtained a predicted structure of SK2/calmodulin and predicted binding sites which were consistent with the literature data. Then we investigated the ligand-protein interaction via virtual mutagenesis, flexible docking, and binding affinity calculation. As a result, the binding energies of mutants have similar trends compared with the EC50 values (R = 0.6 for NS309 in V481 mutants). The results indicate that our protocol can be applied to the drug design of structure unknown proteins. Our study also demonstrates that the integration of in silico approaches is feasible and it facilitates the acceleration of new drug discovery.

Keywords: ligand-protein docking, molecular dynamics simulation, computational drug discovery, SK2, structural prediction, binding site prediction, virtual mutation, pharmaceutical industry

## 1. INTRODUCTION

To discover a new drug is an urgent but time-consuming process (Zhou et al., 2016; Gómez-Bombarelli et al., 2018). In the process of new drug development, in silico approaches have been successfully exploited to perform multiple simulations, such as selecting drug candidates from a database via high throughput virtual screening (Aparna et al., 2014; de Ruyck et al., 2016; Vilar et al., 2016; Imam and Gilani, 2017). The application of in silico approaches not only speeds up new drug discovery, but also collects related information, reveals the mechanisms and proposes new hypotheses. Compared to bench research, computational experiments perform highefficiency simulations which considerably reduce the research time and the cost of experiments (Abel et al., 2017). However, in most cases, in silico approaches played an assisting role in the process. In this study, we propose an all-computational protocol integrating multiple in silico approaches to simulate the entire drug discovery process from de novo protein structure prediction to drug-protein interaction disclosure.

#### Edited by:

Teodorico Castro Ramalho, Universidade Federal de Lavras, Brazil

#### Reviewed by:

Daniel Henriques Soares Leal, Federal University of Itajubá, Brazil Daiana Mancini, Universidade Federal de Lavras, Brazil

\*Correspondence:

Zhong-Ru Xie<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Theoretical and Computational Chemistry, a section of the journal Frontiers in Chemistry

Received: 07 October 2019 Accepted: 24 January 2020 Published: 12 February 2020

#### Citation:

Wu Y, Lou L and Xie Z-R (2020) A Pilot Study of All-Computational Drug Design Protocol–From Structure Prediction to Interaction Analysis. Front. Chem. 8:81. doi: 10.3389/fchem.2020.00081

Up to now, reliable in silico approaches, such as molecular dynamics simulation (MD simulation) and machine learning (Durrant and McCammon, 2011; Lavecchia, 2015; Mortier et al., 2015), have been increasingly developed and applied in finding new drugs and optimizing them for treatment of diseases. However, most research projects only use one single in silico approach. For instance, the homology models of CIB2, a calcium- and integrin- binding protein, were constructed based on CIB1 structures using SWISS-MODEL server (Waterhouse et al., 2018). Based on the models, the way how the point mutations affect the affinities of calcium- and integrin-binding were predicted and then validated by in vitro experiments (Riazuddin et al., 2012). Besides, the results of ligand-protein docking were used to test the substrate specificity of OCT-1 and OCT-2 (organic cation transporter) to guide the following in vivo experimental validation (Papaluca and Ramotar, 2016). These works simply employed the in silico approaches as supportive methods, which did not sufficiently leverage the highspeed advantages of computational methods. Integrating the established in silico researches into an all-computation pipeline and producing validated good results is a milestone in the "omics" era.

Human small conductance calcium-activated (SK2) ion channels, consisted of SK2 subunits and calmodulin molecule, have been proved to be therapeutic targets for treatment of neuronal diseases, such as Parkinson's and amyotrophic lateral sclerosis (ALS) (Bond et al., 2005; Lu et al., 2009). The crystalized structures of human SK2 bound with Riluzole, an approved drug for ALS (Romano et al., 2015), and an analog [pdb (Berman et al., 2003) ligand ID: 658] of the SK2 activator CyPPA, an antiataxic agent (Herrik et al., 2012), have been released recently. It was also reported that ligands of two different chemical classes, Riluzole and its analog NS309, and CyPPA and its analogs, all bind to the same binding site where the interface of SK2 and calmodulin is. In addition, two key residues in the binding site were mutated to investigate how different mutations affect the potency of two ligands (Cho et al., 2018). As a proof of concept, we chose SK2 as the target to examine the protocol we proposed. Hence, in this study, we basically repeated the entire procedure of the previous study of Cho et al. (2018) with consecutive in silico approaches and compared our results with those of the bench experiments. We first constructed the human SK2/calmodulin model (PSK2) using SWISS-MODEL server and docked the ligands Riluzole and CyPPA analog 1 (PDB ligand ID: 658, see section 2.1), onto the predicted binding site. The predicted site is consistent with the reported results (Cho et al., 2018). Then the residues V481 and A477 in the binding pocket were virtually mutated and the mutation effects were assessed via binding energies (MM-GBSA 1GBind) calculation. The results show that the binding energies of mutants have similar trends compared with the EC50 (the concentration of a drug that gives half-maximal response) values (R = 0.6 for NS309 in V481 mutants). Overall, our results suggest that this protocol of in silico approach can provide a systematic prediction on the unknown structures of proteins and potential drugs, and they also demonstrate the ability of in silico approaches to speed up the new drug design process.

## 2. MATERIALS AND METHODS

#### 2.1. Protein Structure Preparation

Having a determined or predicted structure of the drug target protein is the first prerequisite of structure-based drug design. To prove the all-computational protocol is valid, we started our process from structure prediction. To predict the structure of SK2 and calmodulin from Homo sapiens, its amino acid sequence was obtained from the UniProt website (uniprot ID: Q9H2S1). Then we uploaded the sequence onto the SWISS-MODEL server, the most widely-used and reliable structure prediction server, to build a structure model (Bienert et al., 2016; Waterhouse et al., 2018) and selected the 3-D complex structure of Rattus norvegicus SK2 ion channels with NS309 (PDB ID: 4J9Z) (Zhang et al., 2013) as the template of structure prediction. 4J9Z was downloaded from RCSB's Protein Data Bank (Berman et al., 2003), and the predicted structure of SK2 and calmodulin were combined in Maestro (11.5 version, Schrodinger). To test the accuracy of modeling, we uploaded the predicted protein structure and the crystalized complex structure of human SK2 and calmodulin with Riluzole (PDB ID: 5V02) onto Zhang's web server (Zhang and Skolnick, 2005) to calculate the TM-align score. 5V02 was downloaded from RCSB's Protein Data Bank. In addition, sitedirected mutants were constructed using the Mutate Residue function of Maestro.

A PDB structure of the target protein cannot be directly used in molecular docking without preprocessing. In most of the cases, a PDB file does not include the information of hydrogens, the (potential) charges of atoms, or the bond orders between any two atoms. In addition, the protein structures may be determined with a missing fragment(s), a low resolution or alternate positions, or under an unnatural condition, for example, low or high pH values. To make sure molecular docking can simulate the binding between ligands and the target protein correctly and precisely, the protein and ligand preparation is necessary. The wild-type and mutant (predicted) protein structures to be used for docking were processed by protein preparation wizard in Maestro (Sastry et al., 2013). The workflow of protein preparation contains three steps as follows: (1) Preprocess: assigning bond orders, adding hydrogens, creating zero-order bonds to metals, creating disulfide, filling in missing side chains using Prime, deleting waters beyond 5.00 Å from het groups (to keep the water molecules which may form hydrogen-bond bridges between the protein and the ligand and remove those cannot form hydrogen-bond bridges) and generating het states using Epik (pH = 7.0 ± 2.0) (Shelley et al., 2007); (2) Optimization: setting pH = 7.0 and performing optimization; (3) Minimization: this step was performed using the OPLS3 force field (Harder et al., 2015). The converge heavy atoms to root-mean-square deviation (RMSD) is 0.30 Å.

#### 2.2. Ligand Preparation

The 3-D molecular structures of Riluzole and NS309 were obtained from the PubChem database (Kim et al., 2018). The 3-D molecular structures of CyPPA analog 1 ((4-chloro-phenyl)- [2-(3,5-dimethyl-pyrazol-1-yl)-pyrimidin-4-yl]-amine) and analog 2 (4-chloro-phenyl)-[2-(3,5-dimethyl-pyrazol-1-yl)-6 methyl-pyrimidin-4-yl]-amine) were built in Maestro (11.5 version, Schrodinger) based on a previous study (Cho et al., 2018). All the compounds were prepared using OPLS3 force field in Ligprep panel in Maestro (Sastry et al., 2013; Harder et al., 2015). The preparation process included converting 2D structures to 3D ones, adding hydrogens, computing correct partial charges, and optimizing the structures.

### 2.3. Binding Site Prediction

Knowing the potential ligand binding site(s) is also an important prerequisite prior to molecular docking. There are many welldeveloped binding site prediction methods and servers (Xie and Hwang, 2015); however, the predictions produced by different methods may disagree with each other. Therefore, researchers usually compare the prediction results produced by different methods to find the consensus among all predictions. In our study, binding site prediction process was completed using the LISE web server (http://lise.ibms.sinica.edu.tw/applet/) (Xie and Hwang, 2012), which is reported to have the highest accuracy (80–90% for a soluble protein), and the binding site prediction tool, Sitemap, in Maestro, which we used in the docking procedure after this step. The SK2 predicted structure was uploaded onto the LISE web server for binding site prediction. The top three predictions from LISE were then downloaded and imported into Maestro. Meanwhile, in Maestro, SiteMap (Halgren, 2007, 2009) predicted five ligand binding sites. After comparing the results obtained using two prediction methods, we used the consensus region to define the docking grid box.

#### 2.4. Ligand-Protein Docking

In order to predict the details of the interaction between ligands and the target protein and to estimate their binding affinities (see section 2.5), ligand docking was conducted the extra-precision (XP) mode in Maestro. Maestro has three precision options for docking including high throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP). Users can choose an option based on their need or the computational load. After the ligands and the target proteins were processed using Ligprep and protein preparation, respectively, a receptor grid box was generated according to the results of binding site prediction. The size of the receptor grid box was set as default. To investigate the interaction of the protein and ligands, Induced Fit Docking (IFD) (Farid et al., 2006; Sherman et al., 2006a,b; Clark et al., 2016) was performed in Maestro. Using the IFD, the Ligprep outputs were imported and docked to the target protein. The standard protocol was applied to generate up to 20 poses. The force field was OPLS3. Under the prime refinement tab, the conformations of binding site residues within 5 Å (default value) of the ligand were refined. In the Glide redocking process, the energy of each protein/ligand complex structure and the number of top structures were set as the default settings. The XP mode was used for all IFD process.

#### 2.5. MM-GBSA Calculation

The binding energy (1GBind) between a protein and a ligand reflects how stable they bind to each other and how a point mutation affects the ligand binding. Therefore, we examined if our model can correctly predict the trend of binding affinity changes of the mutations on the target protein. In this study, 1GBind were estimated using the Prime MM-GBSA module in Maestro (Greenidge et al., 2012). In MM-GBSA panel, the pose viewer files of docked complex were uploaded. The solvation model was VSGB and the force field was OPLS3 (Li et al., 2011). Prime MM-GBSA 1GBind was calculated using this equation:

$$\begin{aligned} \Delta \mathbf{G}\_{\text{Bind}} &= \mathbf{E}\_{\text{complex}} \text{(minimized)} - \left[ \mathbf{E}\_{\text{ligand}} \text{(minimized)} \right] \\ &+ \mathbf{E}\_{\text{receptor}} \text{(minimized)} \end{aligned} \tag{1}$$

Where 1GBind is binding free energy and Ecomplex(minimized), Eligand(minimized), and Ereceptor(minimized) are minimized energies of receptor-ligand complex, ligand and receptor, respectively.

#### 2.6. Virtual Mutation

Based on the literature, valine 481 (V481) of SK2 is a crucial residue which forms the hydrophobic core between SK2 and calmodulin (Cho et al., 2018). To investigate the impacts of V481 mutations in the binding pocket using in silico approaches, we first implemented the site-directed mutagenesis in Maestro. The V481 was virtually mutated to alanine, serine, threonine, aspartate, or phenylalanine. Then NS309 and CyPPA analog 2 were docked onto the mutated binding site of PSK2 using IFD and the binding free energies were calculated using MM-GBSA to reveal the effects of mutated residues.

Alanine 477 (A477) is another vital residue in the binding pocket (Cho et al., 2018). We exploited the same method mentioned above to virtually mutate A477 to isoleucine, leucine, valine, serine, threonine, arginine, and aspartate, docked NS309 and CyPPA analog 2 onto the mutated binding site of PSK2 and calculated the binding free energies using MM-GBSA.

#### 2.7. Molecular Dynamics Simulation

The molecular dynamics (MD) simulations were performed using GROMACS version 2018.1 and CHARMM36 all-atom force field (March 2019) (Vanommeslaeghe et al., 2010, 2012; Vanommeslaeghe and MacKerell, 2012; Yu et al., 2012; Gutiérrez et al., 2016). The starting coordinates of each protein-ligand complex were obtained from docking experiments. Then we defined a dodecahedral unit cell and filled it with water molecules. After adding ions, the complex was minimized for 50,000 steps of steepest descent minimization. Next, the complex was equilibrated using an NVT ensemble (constant Number of particles, Volume, and Temperature) and NPT ensemble (the Number of particles, Pressure, and Temperature). The target temperature for equilibration was 300 K. At last, the simulations were performed for 30 ns. After the MD simulations, we calculated the RMSD of the residues which were mentioned in the previous research in four trajectories (Cho et al., 2018). Then, we selected four time points of two residues: I100 on protein-NS309 complex (15,000, 18,000, 24,000, and 30,000 ps) and D64 on protein-CyPPA analog 1 (2,610, 6,000, 15,000, and 21,000 ps). Finally, four conformations of both residues were converted into PDB files and were superposed using PyMol.

### 3. RESULTS

### 3.1. Structure Prediction of Human SK2/Calmodulin Complex

The structure of human SK2/calmodulin was generated using template-based modeling on the SWISS-MODEL server (**Figure 1A**). Based on the structure of SK2/calmodulin from Rattus norvegicus (PDB ID: 4J9Z), the structure models were predicted using the amino acid sequences of human SK2 and calmodulin. Additionally, the loop (residue A403 to residue D413, the "intrinsically disordered fragment"—IDF), which had not been determined in structure of human SK2/calmodulin complex (PDB: 5V02), was obtained in the predicted model (**Figure 1A**). It suggests that the predicted structure can be used to supplement the crystallized structure.

To test the accuracy of the predicted model, we used TMalign server to calculate the TM-align score (Zhang and Skolnick, 2005). The predicted models were aligned with 5V02, and the TM-align scores of SK2 and calmodulin are 0.99124 and 0.89215. These TM-align scores show that the structure of human SK2/calmodulin complex has been accurately predicted.

## 3.2. Binding Site Prediction

To determine the binding pocket in SK2/calmodulin complex, the Sitemap and LISE were exploited to predict binding sites (Halgren, 2007, 2009; Xie and Hwang, 2012). The top three predicted results from LISE were overlapped with the results of

FIGURE 1 | Overlapped structures of the predicted model of human SK2 ion channels and 5V02. (A) The predicted structure of human SK2 ion channels (PSK2) in blue overlapped with its crystallized structure (5V02) in red. The structure of a missing loop (IDF) near the ligand binding pocket in 5V02/5V03 has also been predicted. The green compound was Riluzole, which was the ligand in 5V02. (B) The predicted binding sites (Blue dots represent the results from Sitemap; red dots represent the results from LISE) overlapped with the Riluzole binding site in 5V02.

Sitemap. After the comparison, we found that there was only one consensus. Then we selected this binding site to generate receptor grid for docking. To verify the accuracy of binding site prediction, we also overlapped the predicted binding site with that of 5V02 (**Figure 1B**). Notably, the predicted binding site is the same binding site of Riluzole in 5V02, which suggests that this binding site is the potential binding site for the ligands. Hence, the in silico approaches successfully predict the accurate binding site.

## 3.3. Molecular Docking

Based on the previous study (Cho et al., 2018), the interface of SK2 and calmodulin can be bound by Riluzole, NS309, CyPPA analog 1, and analog 2. To investigate whether we can obtain the same results using in silico approaches, we first docked Riluzole and CyPPA analog 1 onto the predicted model (PSK2) via IFD (Induced Fit Docking) in Maestro, and redocked these ligands onto the determined structures (5V02 and 5V03) and calculated the binding energies as the control. Then we calculated the binding energies using MM-GBSA to estimate the binding affinities (Greenidge et al., 2012). As a result, the MM-GBSA 1GBind of PSK2 with Riluzole and CyPPA analog 1 are −40.19 and −56.11 kcal/mol, respectively. As indicated in **Table 1**, those results of PSK2 are compatible to the results of 5V02 and 5V03, which demonstrate that accurate results can be obtained using in silico approaches. Additionally, the docking pose of PSK2 with Riluzole is almost identical with that in 5V02 (**Figure 2A**, the RMSD between the ligand of crystal structure and docking poses on 5V02 or PSK2 is 0.43 or 0.72 Å), which indicates that the simulated result from IFD can obtain accurate data in comparison with the results of crystallization. In **Figure 2B**, the coordinates of docked and native ligands are almost the same even though the poses of two ligands are not completely overlapped (the RMSD between the ligand of crystal structure and docking poses on 5V03 or PSK2 is 1.59 or 0.87 Å). Analyzing the docking results, Riluzole and CyPPA analog 1 interact with residues in both SK2 and calmodulin (**Figure 3**). The interacting residues in the binding site are mostly hydrophobic residues. As CyPPA analogs are larger molecules, they interact with more residues. Compared the list of interacting residues (**Table S1**), our results are almost identical to those of the previous study (Cho et al., 2018). The discrepancy may be because Maestro analyzes the interactions between ligands and the protein based on the interaction energy and the previous

TABLE 1 | MM-GBSA 1GBind of ligands bound to crystallized structures and predicted structures.

study (Cho et al., 2018) simply lists the residues within 5 Å of either ligand. According to the docking results of PSK2 with Riluzole and CyPPA analog 1, we docked NS309 and CyPPA analog 2 on the same binding site using the same methods mentioned above (**Figure 4** and **Figure S1**). The binding affinities MM-GBSA 1GBind are showed in **Table 1**. The MM-GBSA 1GBind values of CyPPA analog 1 and CyPPA analog 2 on the PSK2 are more negative than those of the other ligands, which suggests that CyPPA analog 1 and CyPPA analog 2 are promising drug candidates.

5V03 and the red-colored CyPPA analog 1 is the determined structure in 5v03.

To verify the ligand-induced perturbations, NMR spectrum was used to identify residues with significant chemical shifts in previous experiments (Cho et al., 2018). With computational approaches, we ran MD simulations for each protein-ligand complex to simulate the conformational changes after the binding of ligands. A previous research reports that the residues on EF hands of calmodulin had conformational changes due to the ligand binding (Cho et al., 2018). Therefore, we calculated

CyPPA analog 2 (B).

RMSD for those residues. According to the results of RMSD, we compared four trajectories of each residue and superposed their different poses at different time points. As a result, we found I100, on the complex of NS309, had obvious perturbations on four time points (**Figure 5** and **Figure S2**). In addition, D64 which was a critical residue for calcium sensing, also showed dramatically conformational changing on the complex structure of PSK2 and CyPPA analog 1 (**Figure 6** and **Figure S3**). Hence, the different poses of those residues demonstrated that the ligand could induce the perturbations of the calmodulin, which was consistent with the conclusion in the previous research (Cho et al., 2018).

### 3.4. Virtual Mutations at V481 and A477 in the Binding Pocket

The results of a previous study (Cho et al., 2018) show that the site-directed mutations on the key residues V481 and A477 in the binding site result in changes in the binding affinities. To validate that the in silico approach can simulate and predict the impacts of these mutations, we performed the virtual mutation experiments, docked the ligands on all mutants, and calculated the corresponding binding energies of each mutant. The MM-GBSA 1GBind of V481 mutants are shown in **Table 2**. To verify the accuracy of in silico approaches, the Pearson's correlation coefficients between MM-GBSA 1GBind and EC50 were calculated (R = 0.6 for NS309 in V481 mutants). **Table 2** indicates that substitutions with small side chains (V481S) or a charged amino acid (V481D) significantly decrease the binding affinities of NS309. Conversely, the binding affinity of NS309 in PSK2 V481F is close to that in wild-type. The results above are consistent with the conclusion in literature, that is, this position requires a non-charged residue with a bigger side chain (Cho et al., 2018).

Similarly, we find that CyPPA analog 2 in PSK2 V481F mutant with large aromatic side chains also shows a close binding affinity compared to that in PSK2 WT (**Table 2**), which demonstrates a good correlation between calculated binding affinities and EC50. The mutants with small side chains (V481A and V481S) or a charged amino acid (V481D) also shows the relatively lower binding affinities of CyPPA analog 2. Those results demonstrate that the simulation results are compatible with the data from biological experiments. The consistent conclusion has successfully proved that the all-computational protocol can be widely applied in future biomedical research.

In **Table 3**, all PSK2 A477 mutants have slightly lower NS309 potency than that of PSK2 WT and the predicted binding affinities MM-GBSA 1GBind have similar results. According to the results in literature (Cho et al., 2018), CyPPA analog 2 should exhibit shifted potencies at PSK2 A477L, PSK2 A477V, PSK2 A477S, PSK2 A477T, PSK2 A477R, and PSK2 A477D, but not at PSK2 A477I. However, in **Table 3**, the MM-GBSA 1GBind of PSK2 A477I is not different from those of other mutants. As an isoleucine has a long side chain, different rotamers may largely change the estimated binding affinities. Performing an MD simulation may be a good solution to optimize the structures of mutants and improve the docked poses and estimated binding energies.

## 4. DISCUSSION

In the field of new drug discovery, research efficiency is particularly essential. On one hand, the speed of new drug development is of great importance to patients, especially the ones with fatal diseases such as cancers or acute infectious diseases (Shi et al., 2015). Statistics show that there will be around 1.7 million new cancer cases and 600 thousands cancer deaths in the United States in 2019 (Siegel et al., 2019). This race against time has always been a huge challenge for the researchers in this field. On the other hand, to bring a new drug to the market from compound identification to final FDA approval usually costs up to billions of dollars (Cleary et al., 2018). Therefore, the time and cost-efficient virtual process of drug development will benefit many pharmaceutical companies and our society. In silico approaches which can save considerable amount of research time and cost should be applied to drug design. With the rapid development of computer science and engineering, the availability and accuracy of in silico approaches have been constantly improving. Although many progresses have been made in utilizing in silico approaches to simulate certain biological experiments, the whole experimental processes completed in the virtual way from protein structure simulation

to protein-drug interaction characterization has never been achieved before.

In this study, we proposed a protocol of integration of in silico approaches to simulate the process from protein structure determination to key residues mutagenesis and characterization. To validate this strategy, we selected the human SK2 ion channels as our target protein. With successful prediction, we obtained accurate protein structures (TM-align score >0.5) and the same binding site as the crystallized structure. Furthermore, the docking poses of Riluzole and CyPPA analog 1 are consistent with the ligand bound conformations in the crystallized structures. We also successfully reproduced similar effects of site-directed mutagenesis on the ligand binding, which demonstrated great potential of the integration of in silico approaches. However, the purpose of integrating in silico approaches is not to completely replace biological experiments but to speed up the drug discovery process with the continuous and automatic computational process.

A possible reason why an all-computation protocol of drug design has not been proposed and implemented is the inaccuracy or uncertainty of prediction results might accumulate in the sequential computational pipeline. However, the state-of-the-art bioinformatics algorithms, software or servers have been highly accurate in many cases so that it is time to integrate them into a fully computational process or even a fully automatic process. This is the first study to demonstrate the feasibility of an all-computational protocol in drug design. To achieve the ultimate goal of "automatic" drug discovery, more online servers or computational algorithms like PROCHECK (Laskowski et al., 1993), which can be used to assess or estimate the reliability of prediction results generated by each in silico approach, need to be developed. Despite its innovative approach, there are a few limitations of this study. First, the accuracy of protein structure prediction relies on the methods and/or templates. In our research, we selected the SK2/calmodulin from Rattus norvegicus (PDB ID: 4J9Z) as the template to build protein structure on SWISS MODEL, whose results are more accurate than the results from other webservers (data not shown). Hence, a reliable tool or method is critical to the accuracy of the final simulated results. Second, the binding affinity changes are hard to be precisely reproduced, especially those of the mutants, because considering possible conformational changes on target proteins is still the biggest challenge of docking. This suggests that MD simulation which can simulate the conformational changes should be used to further improve the precision of the predictions.

In conclusion, this work established and demonstrated an integrated protocol of in silico approaches for the first

#### REFERENCES

time. Its applicability can be potentially extended beyond the characterizing of SK2 ion channels to investigating other proteins with unknown structures, such as the Alzheimer's disease related proteins (Fitzpatrick et al., 2017; Hatami et al., 2017), which are also treatment targets of neural degenerative diseases. Although there are challenges to the in silico approaches, our work still paves a new way toward an automatic procedure of drug design.

#### DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/**Supplementary Material**.

### AUTHOR CONTRIBUTIONS

YW carried out the research, analyzed the data, and drafted the manuscript. LL analyzed the data. Z-RX conceived the project, guided the research, and revised the manuscript. All authors read and approved the final manuscript.

### FUNDING

This work has been supported by a start-up grant from the College of Engineering, University of Georgia.

#### ACKNOWLEDGMENTS

The authors would like to acknowledge Dr. Gareth Young and Dr. Shenping Liu for kindly providing their experimental data (EC50), and Dr. Wen-Ching Chan, Xingzi Yuan, Hsin-Tzu Wang, Jane Hua, and Crystal Zhu for the useful suggestions on editing this manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2020.00081/full#supplementary-material

Durrant, J. D., and McCammon, J. A. (2011). Molecular dynamics simulations and drug discovery. BMC Biol. 9:71. doi: 10.1186/1741-7007-9-71

randomised, double-blind, placebo-controlled trial. Lancet Neurol. 14, 985– 991. doi: 10.1016/S1474-4422(15)00201-X

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Wu, Lou and Xie. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# In silico Strategies to Support Fragment-to-Lead Optimization in Drug Discovery

Lauro Ribeiro de Souza Neto1†, José Teófilo Moreira-Filho2†, Bruno Junior Neves 2,3† , Rocío Lucía Beatriz Riveros Maidana1,4, Ana Carolina Ramos Guimarães 4 , Nicholas Furnham5 \*, Carolina Horta Andrade2 \* and Floriano Paes Silva Jr. 1 \*

#### Edited by:

Teodorico Castro Ramalho, Universidade Federal de Lavras, Brazil

#### Reviewed by:

Hyun Lee, University of Illinois at Chicago, United States Salvatore Guccione, University of Catania, Italy Daiana Mancini, Universidade Federal de Lavras, Brazil

#### \*Correspondence:

Nicholas Furnham<EMAIL_ADDRESS>Carolina Horta Andrade<EMAIL_ADDRESS>Floriano Paes Silva Jr<EMAIL_ADDRESS>
†These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Medicinal and Pharmaceutical Chemistry, a section of the journal Frontiers in Chemistry

Received: 25 October 2019 Accepted: 30 January 2020 Published: 18 February 2020

#### Citation:

de Souza Neto LR, Moreira-Filho JT, Neves BJ, Maidana RLBR, Guimarães ACR, Furnham N, Andrade CH and Silva FP Jr (2020) In silico Strategies to Support Fragment-to-Lead Optimization in Drug Discovery. Front. Chem. 8:93. doi: 10.3389/fchem.2020.00093 1 LaBECFar – Laboratório de Bioquímica Experimental e Computacional de Fármacos, Instituto Oswaldo Cruz, Fundação Oswaldo Cruz, Rio de Janeiro, Brazil, 2 LabMol – Laboratory for Molecular Modeling and Drug Design, Faculdade de Farmácia, Universidade Federal de Goiás, Goiânia, Brazil, 3 Laboratory of Cheminformatics, Centro Universitário de Anápolis – UniEVANGÉLICA, Anápolis, Brazil, 4 Laboratório de Genômica Funcional e Bioinformática, Instituto Oswaldo Cruz, Fundação Oswaldo Cruz, Rio de Janeiro, Brazil, 5 Department of Infection Biology, London School of Hygiene and Tropical Medicine, London, United Kingdom

Fragment-based drug (or lead) discovery (FBDD or FBLD) has developed in the last two decades to become a successful key technology in the pharmaceutical industry for early stage drug discovery and development. The FBDD strategy consists of screening low molecular weight compounds against macromolecular targets (usually proteins) of clinical relevance. These small molecular fragments can bind at one or more sites on the target and act as starting points for the development of lead compounds. In developing the fragments attractive features that can translate into compounds with favorable physical, pharmacokinetics and toxicity (ADMET—absorption, distribution, metabolism, excretion, and toxicity) properties can be integrated. Structure-enabled fragment screening campaigns use a combination of screening by a range of biophysical techniques, such as differential scanning fluorimetry, surface plasmon resonance, and thermophoresis, followed by structural characterization of fragment binding using NMR or X-ray crystallography. Structural characterization is also used in subsequent analysis for growing fragments of selected screening hits. The latest iteration of the FBDD workflow employs a high-throughput methodology of massively parallel screening by X-ray crystallography of individually soaked fragments. In this review we will outline the FBDD strategies and explore a variety of in silico approaches to support the follow-up fragment-to-lead optimization of either: growing, linking, and merging. These fragment expansion strategies include hot spot analysis, druggability prediction, SAR (structure-activity relationships) by catalog methods, application of machine learning/deep learning models for virtual screening and several de novo design methods for proposing synthesizable new compounds. Finally, we will highlight recent case studies in fragment-based drug discovery where in silico methods have successfully contributed to the development of lead compounds.

Keywords: fragment-based, drug discovery, lead discovery, in silico methods, machine learning, de novo design, optimization, hot spot analysis

## INTRODUCTION

#### Fragment-Based Drug Discovery

Since the inception of fragment-based drug discovery (FBDD) over 20 years ago it has become an established technology used in both industry and academia (Hubbard, 2015). FBDD offers an attractive approach for effectively exploring the chemical space for binding a target protein. In conventional high-throughput screening (HTS) campaigns, large libraries of often complex compounds are screened for activity against a target (Hall et al., 2014). In contrast, FBDD use relatively small libraries of low complexity compounds representing fragments of larger more drug-like compounds. By reducing the complexity of the chemicals screened more of the potential binding sites of a target protein can be explored through the binding promiscuity of the fragments (Thomas et al., 2017). Where fragments do bind, albeit with lower potency than the drug-like molecules of HTS, they offer good starting points to design larger higher affinity binders using knowledge of the protein structure as a template to generate compounds with greater ligand efficiency (improved per atom binding energy to the target). This bottom-up approach means that a greater range of chemical space can be explored, leading quickly to higher affinity lead compounds with greater specificity (Patel et al., 2014).

FBDD projects require relatively lower investments in research and development (R&D) than HTS (Davis and Roughley, 2017). An example is the discovery of vemurafenib (ZelborafTM), the first fragment-derived drug, which moved relatively very quickly (6 years) between the phases of R&D pipeline before reaching Food and Drug Association (FDA) approval (Erlanson et al., 2016). Thus, FBDD provides attractive opportunities for the drug discovery field.

### Output of Structure-Enabled Fragment Screening Campaigns

FBDD workflows are multi-step starting with target selection and protein isolation and followed by an initial screen of the fragment library using biophysical techniques such as nuclear magnetic resonance (NMR), surface plasmon resonance (SPR), thermal-shift assay, microscale thermophoresis (MST), mass spectrometry, and others. For fragments which show evidence of binding, a further step of hit validation and characterization occurs principally using X-ray crystallography (Verdonk and Hartshorn, 2004). Using hit characterization, an iterative cycle of fragment development can occur employing a range of in silico and experimental techniques. Advances in this protocol try to compress the process by combining the initial fragment screen with the hit characterization. This has been implemented in a high throughput FBDD platform called XChem located at the United Kingdom's national synchrotron the Diamond Light Source (Cox et al., 2016). It uses the ability to produce and handle a large number of crystals of the target protein to screen the fragment library by soaking each individual crystal with a fragment and then using X-ray crystallography to determine which fragments have bound and where. Though this high throughput technique often provides multiple hits, care needs to be taken in interpreting the significance of the hit. Promiscuous fragments may bind parts of the protein which are not involved in the protein function and therefore are unlikely to yield a successful inhibitor. Additionally, as the fragments are by their very nature weak binders and X-ray crystallography being a sensitive technique, observed binding events might be transient and not easily reproducible. It is therefore important to confirm hits with orthogonal structural (e.g., NMR), biophysical techniques (SPR, MST, etc.) or in vitro biological assays.

### Fragment Libraries

A crucial step in FBDD process is in the development and choice of the fragment library used in the screening campaign. Several fragment libraries have been developed that exploit certain properties or chemistries. An example of a fragment library is the Diamond-SGC Poised Library (DSPL) (Cox et al., 2016). This has been developed for use with high-throughput XChem platform and consists of around 760 fragments that have been selected to contain at least one functional group that is open to rapid, cheap follow-up synthesis using robust well-characterized reactions (poised) and maximizing chemical diversity. Other fragment libraries optimize other properties such as solubility, 3D traits or based on subsets of existing drugs and related molecules such as natural products (Schuffenhauer et al., 2005). The fragment libraries generally share similar properties of "Rule of 3" compliant i.e., less than 300 Da molecular weight, 3 or less hydrogen bond donors, 3 or less hydrogen bond acceptors and CLogP no more than 3. In addition, they are soluble in dimethyl sulfoxide (DMSO) or phosphate buffered saline. Fragment libraries generally tend to be <1,000 fragments, which is significantly less than the many millions of compounds screened in high-throughput and high content screening campaigns (Trevizani et al., 2017).

### Fragment Expansion Strategies

Once the fragment screen has been completed and hits characterized, the next step is the challenge of expanding these fragments to generate larger molecular entities with high binding affinity and demonstrating inhibition activity. There are several strategies that can be followed (Lamoree and Hubbard, 2017) (**Figure 1**). One option is to use expert medicinal chemistry advice to design and synthesize larger molecules based on the protein and the fragment pose. Another approach is to define vectors along the fragment molecule based on the steric hindrance of the protein target in which the fragment can be expanded. The fragment is then searched for within a large library of synthesizable (or purchasable) molecules which are bigger by between one and three heavy atoms along the identified vectors. These expanded fragments can be synthesized and soaked/cocrystallized and re-screened by X-ray crystallography. Expanded fragments that show improved binding can be further extended or structurally modified using the same process, with this cycle continuing until a larger high-affinity binding entity is reached.

An alternative to this "small steps" approach is to try to get a larger higher affinity binding molecule in a single step. This can be achieved by having an in silico method using the fragment in a substructure search of a large purchasable compound library (e.g., Zinc15), and to virtually screen the results using the pose

Regardless of the route, expanded fragments should be checked for biological activity using in vitro, ex vivo, or in vivo assays.

of the fragment to dock the molecules and rank them based on docking parameters (Trevizani et al., 2017). The top-ranking virtual screening hits can then be co-crystallized and well as evaluated in vitro and in vivo. A final expansion option is to link or merge fragments that hit near to each other or within the same site (Davis and Roughley, 2017). The combined fragments can then be further expanded using the approaches described previously. It is vital that as expansion progresses in vitro and in vivo assays are conducted to asses activity of the new molecules.

In the next section, we will discuss in depth the main optimization approaches used for a fragment structurally characterized in a binding site of its target. Further sections will describe existing software tools or modeling techniques (e.g., machine learning) employed for taking a fragment hit thorough the path for becoming a lead compound—a process known as fragment-to-lead (F2L)—for drug development and conclude by presenting case studies where in silico strategies have been successfully utilized to support the F2L optimization process.

#### FRAGMENT OPTIMIZATION APPROACHES

After the hit identification in a FBDD campaign, the fragment moves forward to the optimization phase. This optimization takes into account the structural characteristics of the ligand as well as its binding site. The principle in using fragments relies on the premise that these molecular entities are more efficient ligands compared to drug-like molecules, and their structures can be further optimized more efficiently. In fact, this constitutes one of its many advantages. As small entities, molecular fragments can be iteratively optimized to show a better pharmacokinetic profile in the later development stages. Druglike molecules may contain functional groups that contribute poorly to protein binding or, in some cases, can even disrupt the protein-ligand interaction. On the other hand, fragments often form high-quality interactions able to more easily bind to the protein target, translating to a greater number of hits. **Figure 2** depicts schematically this concept.

Another advantage of FBDD is the potential for faster hit progression through the campaign, since the fragments are usually structurally simple and many follow-up compounds can be easily purchased from commercial databases (e.g., MolPort, ZINC15, and ChemBridge) instead of being synthesized. Another important characteristic often used to defend this approach is the high hit rates. In this sense, high hit rates means that the FBDD yields relatively more hits in comparison to the traditional methods such as HTS (Coutard et al., 2014; Mondal et al., 2015). This is due the inversely related nature between molecular complexity and the binding probability (Hann et al., 2001). Other advantages includes the more efficient chemical space sampling (Coutard et al., 2014; Mondal et al., 2015) and the relative low cost to implement the FBDD, as it can be seen from comparing the usual size of the HTS library (thousands of compounds) with fragment libraries (hundreds of compounds) (Macarron et al., 2011).

Assessment of the interactions between the fragment and its binding site should be carefully performed for further identification of synthetically accessible vectors on the ligand. Although x-ray crystallography data is a valuable technique in fragment optimization, is important to keep in mind that the observed structural data only represents a snapshot of the system under investigation. It's been known that the ligand affinity can be affected by the structural protein dynamics without changing the ligand-binding interface (Matias et al., 2000; Seo et al., 2014). This complex dynamic environment (Henzler-Wildman and Kern, 2007; Boehr et al., 2009) can affect small and weak ligands as

fragments. With this in mind, many methods can additionally be used to guide the fragment identification/optimization either providing complementary data (e.g., thermodynamic data) or acting as orthogonal approaches (Ciulli, 2013a). These methods are mostly biophysical (Shuker et al., 1996; Lo et al., 2004; Navratilova and Hopkins, 2010; Pedro and Quinn, 2016) and their use has some advantages such as, direct measurement of the binding, detection of small ligands with low affinity, and not needing any prior information about the protein function (Ciulli, 2013b). Despite the supremacy of biophysical methods, biochemical approaches are increasingly being used (Godemann et al., 2009; Boettcher et al., 2010; Mondal et al., 2015) in FBLD.

pockets with high complementarity with ligand are highlighted in gray.

In addition to orthogonal and complementary methods, the ligand efficiency (LE) or one of its related metrics should ideally be used to keep track of the quality of follow-up ligands as they progress through the iterative optimization cycle. Some of these parameters are described below.

AHeavy Atom Count; BMolecular Weight; CTopological Polar Surface Area;

For the sake of brevity these metrics will not be further discussed and we recommend the references above for a deeper understanding. The structural complexity of the protein makes larger, more complex and less efficient molecules less likely to bind. This is one of the main reasons why fragment libraries often yield more hits when compared to a drug-like molecule commonly used in HTS (Hann et al., 2001). The use of fragments is a bottom-up approach, starting from less complex molecules with greater binding efficiency and ending up with a larger optimized molecule. As already highlighted, there are three main strategies that can be employed to optimize a ligand found bound in its target surface: linking, merging and growing (**Figure 3**). The next sections are dedicated to discussing in more depth each of them.

#### Growing

Fragment growing (**Figure 3A**) is the strategy most commonly employed during FBLD campaigns. As the name suggests, it consists of modifying the fragments to increase their size. Conceptually this approach is identical to the traditional compound modification methods employed in the optimization of hits from HTS campaigns. This modification occurs through the addition of groups.

A recent paper published by Strecker et al. (2019) is an example of how the growing strategy can be used to improve bind affinity. Using computer-aided drug design (CADD) and synthesis, the authors explored small structural modifications in a previously (PDB: 3U0X) identified compound (**1**) (Ki = 800 µM).

These studies showed that a modification of a fragment phenyl moiety to a naphthyl allowed two new simultaneous π-π interactions, a parallel-displaced with Trp300 and an edge-to-face with His233. This minor modification led to a compound (**2**) with a 3-fold improved binding affinity (Ki = 271µM) (**Figure 4**).

gray. The red and gray colors represent the level of complementarity of ligand with the active site. Pockets with low complementarity with ligand are colored in red; pockets with high complementarity with ligand are highlighted in gray.

This example highlights the use of optimal growth vectors to introduce a rigid group, which led to an increased binding affinity. Alternatively, introduction of a moiety with increased number of rotatable bonds could impact negatively due to the entropic penalty—in the affinity. Although this optimization approach can be computationally aided without further structural data, small modifications—as in the case of the hypothetical flexible moiety addition—can led to great changes in binding mode. When growing fragments is the chosen approach, structural data can be decisive to avoid misinterpretation.

#### Linking

Fragment linking (**Figure 3B**) describes the process of joining two non-competitive fragments (i.e., fragments that bind in two different sub-pockets of the binding site) with a chemical linker or spacer. Although conceptually simple, linking fragments is perhaps the most challenging strategy to implement. Although fragment linking is the most attractive approach in terms of rapid improvement of potency, the design of a linker with suitable flexibility while not disturbing the original binding modes of the fragments, makes it one of the most challenging optimization approaches.

As previously discussed, the introduction of flexible moieties affects these compounds properties and an optimal orientation should always be pursued. In fact, varying the degree of rigidity of a linker for the purpose of conformational restriction of the linked product can be used as a strategy for linker optimization, as it can be seen in Chung and colleagues work (Chung et al., 2009). This work shows how a conversion of oxime linkers into monoamine and diamines interferes with the rigidity and its impact on binding.

Although often neglected, the impact on the ADMET properties should also be taken in consideration. In the case of the linker, that usually adds rotatable bonds to the system (Ichihara et al., 2011; De Fusco et al., 2017), this modification can lead to poor PK features, like low permeability (Veber et al., 2002).

#### Merging