Speeding up Early Drug Discovery in Antiviral Research: A Fragment

Apr 24, 2017 - Hepatitis C constitutes an unresolved global health problem. This infectious disease is caused by the hepatotropic hepatitis C virus (H...
0 downloads 9 Views 2MB Size
Subscriber access provided by UB + Fachbibliothek Chemie | (FU-Bibliothekssystem)

Article

Speeding up early drug discovery in antiviral research: A fragmentbased in silico approach for the design of virtual anti-hepatitis C leads Alejandro Speck-Planche, and M. Natalia Dias Soeiro Cordeiro ACS Comb. Sci., Just Accepted Manuscript • Publication Date (Web): 24 Apr 2017 Downloaded from http://pubs.acs.org on April 25, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

ACS Combinatorial Science is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

Speeding up early drug discovery in antiviral research: A fragment-based in silico approach for the design of virtual anti-hepatitis C leads Alejandro Speck-Planche,a* and M. Natália D. S. Cordeiro,a a

LAQV@REQUIMTE/Department of Chemistry and Biochemistry, University of Porto, 4169-007 Porto, Portugal. ABSTRACT Hepatitis C constitutes an unresolved global health problem. This infectious disease is caused by the hepatotropic hepatitis C virus (HCV), and it can lead to the occurrence of life-threatening medical conditions such as cirrhosis and liver cancer. Nowadays, major clinical concerns have arisen due to the appearance of multidrug resistance (MDR), and the side effects especially associated with long-term treatments. In this work, we report the first multitasking model for quantitative structure-biological effect relationships (mtk-QSBER), focused on the simultaneous exploration of anti-HCV activity and in vitro safety profiles related to the absorption, distribution, metabolism, elimination, and toxicity (ADMET). The mtk-QSBER model was created from a dataset formed by 40,158 cases, displaying accuracy higher than 95% in both training and prediction (test) sets. Several molecular fragments were selected, and their quantitative contributions to anti-HCV activity and ADMET profiles were calculated. By combining the analysis of the fragments with positive contributions and the physicochemical meanings of the different molecular descriptors in the mtk-QSBER, six new molecules were designed. These new molecules were predicted to exhibit potent anti-HCV activity and desirable in vitro ADMET properties. In addition, the designed molecules have good druglikeness according to the Lipinski’s rule of five and its variants. Keywords: Anti-HCV activity; ADMET; contribution; design; fragments; quadratic indices; mtk-QSBER, screening. To whom correspondence should be addressed: Alejandro Speck-Planche E-mail: [email protected]; Fax: +351 220402659 (ASP) Introduction Infectious diseases have afflicted mankind for centuries. Among them, hepatitis C constitutes an unresolved global health problem. This illness is caused by the hepatotropic hepatitis C virus (HCV), and it is characterized by a liver inflammation, which can subsequently progress to life-threatening conditions such as cirrhosis,1 acute liver failure,2 and liver cancer.3 Hepatitis due to HCV has a negative impact in worldwide in terms of incidence and mortality. Nowadays, epidemiology studies estimate that three to four million people are newly infected each year by HCV.4 In addition, around 170 million people have been chronically infected and are at risk of developing liver disease including cirrhosis and liver cancer, while 350,000 deaths have occurred each year due to all HCV-related causes.4 Current antiviral therapies have been very efficient in halting the progression of hepatitis C. However, in chronic stages, this disease needs a long-term treatment, which can lead to the emergence of multidrug resistance (MDR),5 and the prevalence of serious side effects caused by antiviral drugs.6 Altogether, the discovery of new, potent and safe anti-HCV agents is a great challenge for the scientific community. Despite the scientific advances achieved in the establishment and consolidation of highly powerful experimental methods such as high-throughput screening (HTS), drug discovery remains a very complex, time-consuming, and expensive process.7 The vast chemical space to be covered in the search for new drugs (1063 for small to medium-sized molecules),8 and the exponential growth of high quality experimental data from medicinal chemistry and drug discovery programs, have paved the way to the application of in silico models, which have become important pillars in drug research.9 Thus, in the context of antiviral research, chemoinformatic tools have been very useful in predicting and/or discovering anti-HCV compounds. Nevertheless, in general terms, works reporting the application of chemoinformatic models in antiviral research suffer from at least one of the following unfavorable drawbacks. Firstly, some of the models have been derived using small databases of structurally related molecules, where only one disease-related protein has been considered.10-11 Secondly, even if the current models are derived from relatively large and heterogeneous datasets of organic compounds, they are aimed at modeling only activity,12-15 thus neglecting the safety profiles related to absorption, distribution, 1

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 23

metabolism, elimination, and toxicity (ADMET). Finally, most of these models have not included information regarding the reliability of the assays under which the molecules have been tested.10-16 With the aim of overcoming the aforementioned handicaps, new advanced chemoinformatic models have been develop in order to integrate multiple kinds of chemical and biological data, which allow them to simultaneously predict different biological effects of chemicals against many biological targets (proteins, microorganisms, cells, laboratory animals, and humans), and taking into account also the reliability of the experimental information.17-21 Consequently, such chemoinformatic models may play an essential role in multi-target drug discovery,22 which seems to superior with respect to its mono-target counterpart (classical drug discovery). Multi-target drug discovery aims to create therapeutic drugs that can simultaneously modulate multiple targets in complex diseases, circumventing drug resistance resulting from single-target mutations or expression changes, and decreasing the risk of side effects due to a reduced number of drug-drug interactions.22-24 In the field of antiviral research, some of the aforementioned in silico models have combined in vitro biological activity against the human immunodeficiency virus (HIV) with disease incidence in USA counties,25-27 while a more recent work have been focused on the fragment-based drug design and virtual screening of compounds with potentially high anti-HIV activity and desirable ADMET profiles.28 In any case, to the best of our knowledge, there is yet no model capable of jointly predicting anti-HCV activity and ADMET properties. In this work, we introduce the first multitasking model for quantitative structurebiological effect relationships (mtk-QSBER), devoted to performing simultaneous exploration of the antiHCV activity and ADMET properties at the in vitro level. Particularly, the mtk-QSBER model derived here is applied to the fragment-based design and screening of new molecular entities virtually displaying potent anti-HCV activity and desirable ADMET profiles. Materials and Methods Dataset and molecular descriptors The methodology detailing the generation of mtk-QSBER models have been largely reported,21,29 so we will limit ourselves to give a brief outline, highlighting only the most important aspects. All the experimental data were retrieved from CHEMBL,30-32 and the dataset was formed by 29,863 molecules. Each of them was tested by considering at least 1 out of 12 measures of biological effect (me), against at least 1 out of 72 specific targets (bt). At the same time, each experiment contained 1 out of 3 labels of assay information (ai) indicating whether the assay was focused on the study of binding phenomena, functional/physiological responses, or ADMET profiles. In addition, the assays were linked to at least 1 out of 7 different target mappings (tm) specifying the general types of targets against which the assays were performed. Finally, each experiment was associated with a probabilistic factor (pc) denoting the degree of reliability of the assay (experimental condition). In so doing, pc values equal to 0.55, 0.75, and 1 were attributed to the labels of assay reliability named autocuration, intermediate, and expert, respectively. Each combination of the elements me, bt, ai, tm, and pc, defines a specific experimental condition, which can be represented as an ontology of the form cj→(me, bt, ai, tm, pc).33 It should be emphasized that the molecules present in our dataset have not been assayed against all the experimental conditions. At the same time, one molecule could be assayed against more than one experimental condition. Consequently, our dataset turned out to be formed by 40,158 statistical cases of molecules. Each case in the dataset was annotated to belong to 1 of 2 possible classes/categories, namely positive [BEi(cj) = 1] or negative [BEi(cj) = −1]. All the assignments were realized according to certain cutoff values that are depicted in Table 1. Notice that BEi(cj) is a binary variable that embodies the biological effect of the ith molecule under the experimental condition cj.

2

ACS Paragon Plus Environment

Page 3 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

Table 1. Cutoff values used to annotate the molecules as positive Measure of effect (Units) Biological Profile Definition Concentration required for inhibition EC50 (nM) Anti-hepatitis activity of viral replication by 50% in infected cells. Permeability (nm/s) ADMET (absorption) Permeability in cells. PPB (%)

ADMET (distribution) Plasma protein binding Concentration that inhibits 50% of the ADMET (metabolism) enzymatic activity of a cytochrome. Concentration that inhibits 50% of the ADMET (metabolism) enzymatic activity of a transferase. Inhibitory concentration at 50% ADMET (metabolism) against a transporter protein. Inhibition constant of a chemical ADMET (metabolism) against a transporter protein. ADMET (elimination) Intrinsic clearance in hepatocytes.

IC50 (nM)cyp IC50 (nM)trf IC50 (nM)trp Ki (nM)trp CL_int (µL/min/106 cells) t1/2_lm (hr)

Cutoff ≤350 ≥120 ≤75 ≥5000 ≥8000 ≥8000 ≥7800 ≤6

ADMET (elimination) Half life in liver microsomes. ≥0.5 Inhibitory concentration at 50% ADMET (toxicity) ≥7943.28 against an ion channel. Inhibition constant of a chemical ADMET (toxicity) ≥8500 against an ion channel. Cytotoxic concentration causing ADMET (toxicity) ≥40000 death to 50% of viable cells.

IC50 (nM)ich Ki (nM)ich CC50 (nM)

The SMILES codes for all the 40,158 cases were stored in a file of type *.txt, which was manually changed to *.smi. Afterward, the software Standardizer was employed to transform the *.smi file to *.sdf.34 Subsequently, from the *.sdf file, the program QUBILs-MAS35 was employed to calculate the topological descriptors known as bond-based quadratic indices Bq(x).36 These molecular descriptors are based on the application of discrete mathematics and linear algebra theory to chemistry, and they have been widely used in different areas of research in drug discovery.37-39 These molecular descriptors can be calculated according to the following mathematical expression: n

n

Bqk ( x ) = [ X T ][E]k [ X ] = ∑∑ k eab • xa • xb (1) a =1 b =1

In Eq. 1, keab represents the elements of the kth power of the bond adjacency matrix E of the molecular graph. In addition, n is the number of bonds (without considering bond multiplicity), while xa and xb are the physicochemical properties of the bonds a and b, respectively. Such bonds are adjacent between each other, i.e., they are incident on the same atom. The x value that characterizes each bond can be calculated as:

x=

Pf

δf

+

Pg

δg

(2)

In Eq. 2, Pf and Pg are the physicochemical properties of the atoms f and g that form the bond, while δf and δg denote their corresponding vertex degrees. In this context, the vertex degree δ of an atom is expressed as the number of adjacent atoms to which the atom under analysis is covalently bound. Setup of the mtk-QSBER model By inspecting Eq. 1, one can see that Bq(x) indices only take into account the chemical structures of the molecules, and therefore, they are not able to characterize the different aspects (me, bt, ai, tm, and pc) related to the experimental conditions under which the molecules are assayed. But this issue can be easily solved by resorting to the Box-Jenkins approach, just as it has been applied in diverse fields of research within drug discovery.17-21,40 This approach is based on moving averages,41 and the mathematical formalism can be written in the following form:

avg_Bq k ( x)c j =

1 n( c j )

n(c j )

∑ Bq

k

( x)i

(3)

i=1

3

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 23

In Eq. 3, n(cj) is the number of molecules annotated as positive that were assayed by considering the same element of the experimental condition (ontology) cj. For instance, in the case of the element bt, n(cj) is the number of molecules considered as positive cases [BEi(cj) = 1] that have been experimentally tested against the same biological target. The procedure and reasoning embodied in Eq. 3 were applied to the elements me, bt, ai, and tm. Also, in the same equation, Bq(x)i is the bond-based quadratic index of the ith molecule, while avg_Bq(x)cj represents the arithmetic mean of the Bq(x)i indices. After calculating the avg_Bq(x)cj values, a subsequent expression can be used to calculate a new type of bond-based quadratic index for each molecule, that is:

DBq k ( x ) i c j = pc • [ Bq k ( x) i - avg_Bq k ( x)c j ]

(4)

In Eq. 4, DBq(x)icj is a deviation term and constitutes an adaptation of the Box-Jenkins operators. This improved bond-based quadratic index considers both the chemical structure and a specific element of the experimental condition (ontology) cj under which a molecule has been assayed. Only the DBq(x)icj indices (128 in total) were used as molecular descriptors to seek for the best mtk-QSBER model. The dataset formed by the 40,158 cases was randomly divided into two series: training and prediction (test) sets. The training set was employed in the search for the best model, and it contained 30155 cases, 12,655 annotated as positive, and 17,500 as negative. On the other hand, the prediction set was used to assess the predictive power of the model. This set was formed by 10,003 cases, 4183 assigned as positive, and 5820 negative. As to the data modeling technique for setting up the mtk-QSBER model, we employed linear discriminant analysis (LDA), which was based on a forward step-wise procedure.41 The LDA modeling was carried out using the program STATISTICA.42 In this case, for the application of the aforementioned forward step-wise procedure in selecting each variable, we used the default value of the F statistic provided by the program STATISTICA.41 Thus, only the molecular descriptors with F ≥ 1 and a p-value ≤ 0.05 were selected to enter into the model. In addition, a third criterion was applied; any molecular descriptor was considered to be significant if, with its incorporation in the model, the sensitivity and specificity were increased by at least 2% with respect to the variables previously entered. Other default parameters used by STATISTICA were: maximum number of steps equal to 100, effect of force equal to zero, sweep delta (detection of redundancy in the design matrix and evaluation of the estimability of hypotheses) equal to 10-7, and inverse delta (checking for matrix singularity in matrix inversion calculations) equal to 10-12. The general discriminant equation of the mtk-QSBER model can be expressed in the following form: z

SBEi (c j ) = a0 + ∑ bm • [ DBq k ( x) i c j ]m

(5)

m=1

In Eq. 5, a0 is the constant term, while bm represents the coefficients of the molecular descriptors. Here, it should be pointed out that during the setup of the model, the program STATISTICA takes the initial categorical values of BEi(cj) and transforms them into continues score values of biological effects SBEi(cj). Then, after the application of default procedures used by this software with the aim of finding the best discriminant function, each predicted SBEi(cj) score is converted to its corresponding categorical value [Pred-BEi(cj)]. Collinearity between molecular descriptors was also analyzed. In this sense, it should be pointed out that any two variables are collinear (or if they are highly correlated), then, coefficients accompanying the variables can become unstable.43 This means that any calculation result derived from the use of the model can be masked by the effect of having highly correlated variables. Consequently, in the presence of highly correlated variables, a model may be prone to predict noise (random error) instead of the underlying relationship between a set of independent variables (molecular descriptors) and the response variable of biological effect under study [in our case, BEi(cj)]. Usually, collinearity between any two variables is estimated by means of the Pearson’s correlation coefficient (PCC).44 There is no unified consensus regarding the level of collinearity to be avoided. However, PCC > 0.9 (absolute value) has been suggested as an indication of strong collinearity between variables.43,45 We would like to remark that if collinearity occurs, then, it is advisable to perform the orthogonalization of the variables, which can be achieved by applying principal component analysis (PCA),46-48 or the Randic’s orthogonalization approach.49 Taking into consideration all these ideas, we used a cutoff interval -0.9 < PCC < 0.9 as a criterion for lack of collinearity or strong correlation between the molecular descriptors that entered in the mtk-QSBER model. 4

ACS Paragon Plus Environment

Page 5 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

The quality of the model was assessed through the analysis of several statistical indices such as Wilks’s lambda (λ), chi-square (χ2), p-value, which were calculated only for the training set. Additionally, other indices such as sensitivity (percentages of correct classification for positive cases), specificity (percentages of correct classification for negative cases), accuracy (overall percentage of correct classification), Matthews’s correlation coefficient (MCC), and the area under the receiver operating characteristic (AUROC) curves were calculated.19 These last five indices were calculated for both training and prediction sets, serving to determine both the internal performance and the predictive power of the mtk-QSBER model, respectively. The general steps followed to derive the mtk-QSBER model appear illustrated in Figure 1. Results and Discussion Mtk-QSBER model Following the strategy outlined before, the best model found by us is shown below along with the LDA statistical indices:

SBEi (c j ) = 3.54 •10-3 DBq3 ( AR)me - 0.23DBq1 ( AW )bt -1.24 •10-4 DBq7 ( HYD)ai + 2.69 •10-3 DBq1 ( PSA)ai + 9.12 •10-5 DBq6 ( AW )tm + 2.14 N = 30,155 λ = 0.45 χ 2 = 24,044.38 p - value 71% for both positive and negative cases, which confirms the great performance of the mtk-QSBER model, demonstrating the ability of the molecular descriptors to successfully classify/predict most of the data. Specific details of all the metrics of the type %CCC(cj) can be found in the Supplementary Information 3. Another interesting statistical index is the MCC which had a value of 0.91 and 0.9 for training and prediction sets, respectively. These values express the very strong correlation between the observed and predicted values of the categorical variable of biological effect BEi(cj). The AUROC curves were also used in order to analyze the performance of the mtk-QSBER model (Figure 2), and a value of 0.992 was obtained for both training and prediction sets. This value means that our model does not behave as a random classifier, for which the AUROC curve is equal to 0.5. Altogether, the analysis of the diverse statistics reveals the high quality and predictive power of the mtk-QSBER model. Each molecular descriptor in the mtk-QSBER model was examined in terms of its potential collinearity with respect to the other descriptors (Table 2). Notice that the highest PCC value in the correlation matrix is 0.76. Such correlation value is within the cutoff interval -0.9 < PCC < 0.9 explained above. Therefore, we can conclude that there are no strong correlations between the molecular descriptors. 5

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

Table 2. Correlation matrix obtained for the molecular descriptors in the mtk-QSBER model. DBq3(AR)m e DBq1(AW)bt DBq7(HYD)ai DBq1(PSA)ai DBq6(AW)tm Descriptor DBq3(AR)m e

1

0.61

0.76

0.46

0.65

DBq1(AW)bt

0.61

1

0.3

0.48

0.62

DBq7(HYD)ai

0.76

0.3

1

0.13

0.56

DBq1(PSA)ai

0.46

0.48

0.13

1

0.36

DBq6(AW)tm

0.65

0.62

0.56

0.36

1

Interpretation of the molecular descriptors from a physicochemical point of view Most of the in silico models reported in the scientific literature are used for the sole purpose of predicting diverse biological effects (activity, toxicity, etc.) of molecules. No attention is paid to the fact that important physicochemical information can be gathered from these models, giving them a more mechanistic and phenomenological meaning. A correct and insightful interpretation of the molecular descriptors in a model can point medicinal chemists and drug designers in the right direction regarding the physicochemical and/or structural characteristics that a molecule should have in order to improve a defined biological effect (increment of activity, diminution of toxicity, etc.). To accomplish the task of interpreting the molecular descriptors that entered in the mtk-QSBER model, we will rely on their relative influences, estimated by the calculation of the absolute values of the standardized coefficients (Figure 3). The higher the value of the standardized coefficient of a molecular descriptor, the greater will be its influence on the biological effects. In addition, from Eq. 6, we will extract physicochemical information regarding how the different molecular descriptors should be varied in order to enhance the anti-HCV activity and improve the safety profiles (ADMET) of any molecule. The molecular descriptors of type DBqk(x)icj are weighted by different physicochemical properties, and the letter k represents their order, i.e., the topological distance between any two bonds (pairs of atoms) atoms. As mentioned before, these descriptors simultaneously account for the chemical structures of the molecules and specific elements of the experimental conditions under which the molecules were tested. Moreover, notice that the DBqk(x)icj indices embody the probabilistic factor pc, which indicates the degree of reproducibility of the experimental conditions (assays). By inspecting Figure 3, one can see DBq3(AR)m e is the second most influential descriptor in the mtkQSBER model, expressing the increase of the molar refractivity (polarizability) in regions where any two pairs of atoms are placed at a topological distance equal to 3. This molecular descriptor depends on the chemical structure and the measures of the biological effects. In convergence with DBq3(AR)m e, the descriptor DBq6(AW)tm embodies the chemical structure and the target mapping, indicating the augmentation of the atomic weight in regions where any two pairs of atoms are placed at a topological distance equal to 6. The descriptor DBq6(AW)tm has the lowest importance in the mtk-QSBER model. Notice that DBq3(AR)m e and DBq6(AW)tm have a steric nature, and they are constrained in some way by DBq1(AW)bt, which has the highest influence among all the molecular descriptors present in the mtkQSBER model. Thus, DBq1(AW)bt describes the diminution of the molecular weight in regions where any two pairs of atoms are placed at a topological distance equal to 1; DBq1(AW)bt depends on the molecular structure and the biological targets against which the molecules were assayed. Hydrophobicity usually has an impact on the biological profile of any molecule. In our mtk-QSBER model, the descriptor DBq7(HYD)ai characterizes this property, expressing the decrease of the hydrophobic character of a molecule in regions where any two pairs of atoms are placed at a topological distance equal to 7. Thus, DBq7(HYD)ai is the third most important descriptor, and it depends on the chemical structure and the assay information. Finally, in contrast with the diminution of the hydrophobic character described by DBq7(HYD)ai, it should be expected an augmentation of the polar surface area, particularly in regions where any two pairs of atoms are placed at a topological distance equal to 1. Such localized increase in the polar surface area is embodied by the descriptor DBq1(PSA)ai, (the fourth most influential), which depends on the chemical structure and the assay information. Altogether, from a more chemical point of view, the molecular descriptors present in the mtk-QSBER model express that any molecule will enhance its biological effect (increment of the activity and/or improvement of any ADMET property) if both hydrophobic and polar regions are present in the chemical structure. Aromatic and heteroaromatic rings are preferred over their aliphatic counterparts. Atoms with 6

ACS Paragon Plus Environment

Page 7 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

low atomic weights (carbon, nitrogen, oxygen, and fluorine) should prevail in the chemical structure with respect to heavy atoms (chlorine, bromine, sulfur, and phosphorus) although the presence of the latter is advised if not more than two heavy atoms appear in the molecular structure. The limitation in the number of heavy atoms has direct implications in the future druglikeness of any molecule. At the same time, atoms such as oxygen and nitrogen, which are involved in polar interactions, should be dispersed along the chemical structure of the molecules. Quantitative contributions of the fragments to multiple biological effects In the field of chemoinformatics, a topological index is a numerical parameter derived from a molecular graph (abstract representation of a molecule where atoms are considered as vertices and bonds are considered as edges).50 More than 20 years ago, it was demonstrated that any graph-theoretical invariant (topological index) can be expressed as a linear combination of the number of times in which certain substructures/fragments appear in a molecule.51 This envisages the fact that any regression or classification model can be represented as a linear equation involving fragment descriptors.52 Consequently, if one calculates the topological indices of a fragment and substitute them in a defined regression or classification model, the quantitative contribution of that fragment to a biological effect under study can be obtained.21,29 Several previous studies have shown that quantitative contributions can serve as a guide to medicinal and pharmaceutical chemists, enabling the detection of 2D pharmacophores and toxicophores, among other substructural features.21,29,53 Indeed, such fragment contributions can lead to the generation of more accurate and meaningful structure-activity relationships (SAR), and in a wider context, structure-biological effects relationships. It is essential to point out that fragment contributions are not absolute; they will always be restricted to the dataset (chemico-biological space) from which the in silico model was generated. Nevertheless, they will express the general likeness or tendency of the fragments having an influence (positive or negative) in the biological effects of the molecules. The present mtk-QSBER model based on LDA has been developed from the topological descriptors named bond-based quadratic indices, and we have used this model as a tool to determine the relative influence of several molecular fragments on the biological effects by considering 230 experimental conditions, i.e., combinations of the elements me, bt, ai, and tm, with pc having a constant value of 1. Towards that goal, we have thus selected 40 molecular fragments that can be found in many of the molecules of our dataset (Figure 4). For each molecular fragment, the molecular descriptors were calculated by using Eqs. 1, 3, and 4. Afterward, the descriptors of the type DBqk(x)icj for each fragment were substituted in Eq. 6, yielding a total of 230 × 40 = 9200 SBEi(cj) scores. Subsequently, in order to eliminate (or at least to diminish) any possible dependence on the size, all the SBEi(cj) scores were divided by the total number of bonds (without considering bond multiplicity and including hydrogen atoms) of their corresponding fragments. As a result of this operation, the normalized scores [Norm-SBEi(cj)] were obtained. Then, the mean and the standard deviation of all these Norm-SBEi(cj) scores were calculated. Finally, a standardized score was computed according to the following formalism: the mean was subtracted from each Norm-SBEi(cj) score, and the result of this subtraction was divided by the standard deviation. The standardized scores constitute the quantitative contributions of the molecular fragments to the different biological effects. Here, we should recall that according to the results in Table 2, the molecular descriptors present in the mtk-QSBER model are not highly correlated, and therefore, the effect of collinearity is not expected to mask the real value of the quantitative contributions of the fragments to the different biological effects. A sample of the experimental conditions and the corresponding fragment quantitative contributions appear depicted in Tables 3 and 4, respectively. At the same time, all the quantitative contributions calculated for all the biological effects can be found in the Supplementary Information 4.

7

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 23

Table 3. Different experimental conditions used for the calculation of fragment contributions. cja,b

me

bt

ai

tm

c1

EC50 (nM)

Hepatitis C virus subtype 1b

Functional

Non-molecular

c2

EC50 (nM)

Hepatitis C virus (isolate Con1)

Functional

Non-molecular

c3

Permeability (nm/s)

Caco-2 (Homo sapiens)

ADMET

Mammalian cell

c4

PPB (%)

Plasma (Homo sapiens)

ADMET

Non-molecular

c5

IC50 (nM)cyp

Cytochrome P450 19A1

Binding

Protein

c6

IC50 (nM)cyp

Cytochrome P450 3A4

Binding

Protein

c7

IC50 (nM)trf

Fatty acid synthase

Binding

Protein

c8

IC50 (nM)trp

Serotonin transporter

Binding

Protein

c9

Ki (nM)trp

P-glycoprotein 1

Binding

Protein

c10

CL_int (µl/min/10-6cells)

Hepatocytes (Homo sapiens)

ADMET

Mammalian cell

c11

IC50 (nM)ich

HERG

Binding

Protein

c12

CC50 (nM)

L6 cells (Rattus norvegicus)

ADMET

Mammalian cell

c13

CC50 (nM)

HEK293 cells (Homo sapiens)

ADMET

Mammalian cell

a

The symbol cj is used to represent an ontology/condition as result of the combinations of the elements me, bt, ai, and tm, all of them being affected by the probabilistic factor pc. a In all the experimental conditions depicted in this table, pc = 1, which corresponds to the highest degree of reliability of the experimental conditions. The analysis of Table 4 indicates that there is no fragment with positive contribution against all the experimental conditions. Nevertheless, fragments such as F1, F6-F8, F11, F14, F17, F19, F21, F26-F29, F33, and F39 have positive contributions in most of the experimental conditions, making thus them suitable substructure to be considered in the future design of new molecules with high anti-HCV activity and desirable in vitro ADMET properties. In contrast to the aforementioned fragments, others such as F3, F4, F16, F18, F20, F30, F32, F35, and F40 mostly exhibit negative contributions, and therefore, they are detrimental to almost all the biological effects. The remaining fragments have low positive and/or mixed contributions, meaning that they should be analyzed with caution when designing new molecular entities. Here, it should be pointed out that it is very common that positive fragments can be present in molecules annotated as negative. The opposite is also common, i.e., negative fragments in molecules annotated as positive. This is not surprising, and what it means is that the presence of a defined fragment is not a sufficient condition to consider a molecule to belong to one category or another. Only the presence of different fragments and the way in which they are connected to each other will determine whether a molecule will be positive or negative with respect to one or more biological effects. Table 4. Quantitative contributions of the fragments to diverse biological effects. IDa

c1

c2

c3

F1

0.37

0.47

0.51

F2

0.23

0.30

0.34

c4

c5

c6

c7

c8

c9

c10

c11

c12

c13

-0.33 -0.45

1.00

0.44

-0.52

1.16

0.82

0.31

-0.04

0.52

-0.29 -0.38

0.70

0.28

-0.44

0.83

0.57

0.19

-0.07

0.34

F3

-0.21 -0.14 -0.10 -0.77 -0.87

0.29

-0.16 -0.93

0.42

0.14

-0.26 -0.54 -0.10

F4

-0.20 -0.04

0.03

-1.41 -1.62

0.87

-0.10 -1.74

1.15

0.56

-0.30 -0.91

0.04

F5

0.01

0.08

0.11

-0.48 -0.57

0.45

0.06

-0.62

0.57

0.33

-0.03 -0.28

0.12

F6

0.48

0.57

0.62

-0.22 -0.34

1.11

0.54

-0.41

1.27

0.92

0.42

0.07

0.63

F7

0.36

0.45

0.49

-0.34 -0.47

0.98

0.42

-0.54

1.15

0.80

0.30

-0.05

0.50

F8

0.88

1.01

1.06

-0.05 -0.22

1.72

0.97

-0.31

1.94

1.47

0.81

0.34

1.08

F9

0.03

0.11

0.15

-0.62 -0.74

0.60

0.08

-0.80

0.76

0.44

-0.03 -0.35

0.16 8

ACS Paragon Plus Environment

Page 9 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

F10 0.00

0.07

0.10

-0.52 -0.62

0.47

0.05

-0.67

0.59

0.33

-0.04 -0.31

0.11

F11 0.65

0.77

0.82

-0.19 -0.34

1.40

0.73

-0.42

1.60

1.19

0.59

0.16

0.83

F12 0.34

0.44

0.48

-0.36 -0.48

0.97

0.41

-0.55

1.13

0.79

0.28

-0.07

0.49

F13 0.24

0.32

0.35

-0.36 -0.47

0.77

0.29

-0.53

0.91

0.62

0.19

-0.12

0.36

F14 0.96

1.08

1.14

0.02

-0.14

1.79

1.04

-0.23

2.01

1.55

0.88

0.41

1.15

F15 0.05

0.11

0.14

-0.41 -0.50

0.47

0.09

-0.54

0.58

0.35

0.01

-0.22

0.15

F16 -0.47 -0.42 -0.39 -0.89 -0.97 -0.10 -0.44 -1.01 F17 0.86

0.00

-0.21

-0.51 -0.72 -0.39

0.04

1.49

0.93

-0.03

1.65

1.31

0.80

F18 -0.25 -0.18 -0.15 -0.77 -0.87

0.22

-0.20 -0.92

0.35

0.09

-0.29 -0.55 -0.14

F19 0.88

0.96 1.01

1.00 1.06

0.16

0.45 0.33

1.01

-0.05 -0.22

1.71

0.97

-0.31

1.93

1.47

0.80

F20 -0.19 -0.12 -0.09 -0.68 -0.77

0.26

-0.14 -0.82

0.37

0.13

-0.23 -0.47 -0.08

1.07

F21 1.31

1.47

1.54

0.11

-0.11

2.38

1.41

-0.23

2.66

2.07

1.21

0.60

1.56

F22 0.14

0.23

0.27

-0.56 -0.69

0.76

0.20

-0.75

0.93

0.58

0.08

-0.27

0.28

F23 -0.05

0.03

0.06

-0.65 -0.76

0.48

0.00

-0.82

0.62

0.33

-0.10 -0.41

0.07

F24 -0.10 -0.02

0.02

-0.70 -0.80

0.44

-0.04 -0.86

0.58

0.28

-0.14 -0.45

0.03

F25 0.32

0.40

0.43

-0.28 -0.39

0.85

0.37

-0.45

0.99

0.70

0.27

-0.04

0.44

F26 0.42

0.53

0.58

-0.42 -0.57

1.17

0.49

-0.66

1.36

0.95

0.35

-0.08

0.59

F27 0.35

0.47

0.52

-0.49 -0.64

1.10

0.43

-0.72

1.30

0.89

0.29

-0.14

0.53

F28 0.78

0.91

0.96

-0.15 -0.32

1.62

0.87

-0.41

1.84

1.37

0.71

0.24

0.98

F29 0.78

0.90

0.96

-0.16 -0.32

1.61

0.86

-0.41

1.83

1.37

0.70

0.23

0.97

F30 -0.34 -0.27 -0.24 -0.86 -0.96

0.13

-0.29 -1.01

0.25

0.00

-0.38 -0.65 -0.23

F31 0.20

-0.56 -0.70

0.88

0.27

-0.77

1.06

0.69

0.14

F32 -0.22 -0.15 -0.12 -0.75 -0.84

0.25

-0.18 -0.89

0.37

0.11

-0.27 -0.53 -0.11

F33 0.48

0.58

0.62

-0.29 -0.42

1.16

0.55

-0.50

1.34

0.96

0.41

0.03

0.64

F34 0.22

0.33

0.37

-0.54 -0.68

0.91

0.29

-0.75

1.09

0.71

0.16

-0.22

0.38

F35 -0.39 -0.33 -0.30 -0.86 -0.94

0.02

-0.35 -0.99

0.13

-0.10

-0.43 -0.67 -0.30

F36 0.23

0.32

0.35

-0.42 -0.53

0.81

0.29

-0.60

0.96

0.64

0.18

-0.15

0.36

F37 0.08

0.16

0.19

-0.48 -0.58

0.58

0.13

-0.63

0.71

0.44

0.04

-0.25

0.20

F38 0.19

0.26

0.29

-0.33 -0.42

0.66

0.24

-0.48

0.79

0.53

0.15

-0.11

0.30

F39 1.27

1.41

1.47

0.22

2.21

1.36

-0.07

2.45

1.93

1.18

0.65

1.49

F40 -0.48 -0.43 -0.40 -0.88 -0.95 -0.12 -0.44 -0.99 -0.03 Referred to the same experimental conditions represented in Table 3.

-0.23

-0.51 -0.71 -0.40

0.30

0.35

0.03

-0.25

0.36

a

In silico design and screening of potentially efficient and safe anti-HCV molecules Nowadays, the use of in silico models generated from ligand-based approaches and using large and heterogeneous datasets of molecules is mostly restricted to the virtual screening of external chemical libraries. This means that the current computational models are usually treated as “black boxes”, yielding results of predictions without a real phenomenological meaning. Consequently, there has been an imperious need in assessing the reliability of the predictions in an attempt to partially solve this issue. In this context, concepts such as applicability domain have been developed as an alternative to creating more reliable models.54-57 Nevertheless, there are aspects that will be detrimental to the reliability of the predictions performed by any model regardless of the use of any computational methodology, data analysis method, and applicability domain approach. First, the experimental errors associated with the assay data play an essential role because when a model is constructed from a heterogeneous dataset, the biological and chemical data from many research laboratories are often used. Consequently, it is not possible to control each element that influences the uncertainty of the experimental data, which will definitely affect the future predictions. Second, despite 9

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the huge diversity of molecular descriptors that can be employed in creating in silico models, there are no universal descriptors. As a result, molecular descriptors will never be sufficiently useful for all the datasets because they characterize a reduced portion of the chemical diversity and complexity embodied in the molecules. Even if one combines many families of molecular descriptors into a single model, it will be extremely difficult to reliably predict all the molecules. Third, Todeschini and coworkers demonstrated that the more rigorous the applicability domain approach, the smaller the chemicobiological space that a model can cover.57 The opposite assumption is also valid and envisages the fact that is up to the analyst to select a method to define the applicability domain. Fourth, in most of the computational models reported in the literature, each new molecule is predicted only one time for a particular endpoint, and that puts a lot of pressure in defining the reliability of a single prediction through the use of the applicability domain, when according to the previously explained aspects, the prediction may be wrong when compared with any future experimental validation. The fifth and one of the most important aspects is that when an in silico model is used to perform virtual screening, one should keep in mind that the model performs predictions according to predefined values, which correspond to descriptors of the molecules to be predicted. Notice that most of the molecules to be screened may not follow the rules of the model. With this, we are not referring to the applicability domain (in which the predicted molecules may fall or not) but to the fact that in any model, the molecular descriptors have different degrees of priority/importance. Such degree of hierarchy embodies a realistic phenomenological meaning regarding the physicochemical and structural features that need to be primarily varied in order to enhance a biological effect under study. Consequently, those molecules whose molecular descriptors comply with such degree of hierarchy (with respect to the corresponding descriptors of the compounds in the training set) may have a better chance of being correctly predicted, and the theoretical results should be expected to converge to a greater with the corresponding experimental validation. This fifth aspect is very often neglected regardless of the applicability domain algorithm that is applied to a model. Bearing in mind all these previous ideas, here, rather than using our mtk-QSBER model to perform virtual screening of external datasets, we apply it to rational lead discovery, demonstrating that at least from a computational point of view, new molecular entities virtually exhibiting potent anti-HCV activity and desirable in vitro ADMET properties can be designed. In doing so, we first inspected the molecular fragments and their respective quantitative contributions reported in Figure 4 and Table 4, respectively. We chose only those molecular fragments with positive contributions against all (or almost all) the biological effects as reported in Table 4 and the Supplementary Information 4. The selected molecular fragments were connected or fused to each other in such a way that the resulting designed molecules (Figure 5) obeyed the physicochemical and structural information provided by the different molecular descriptors in the mtk-QSBER model. Six new molecules were designed. Subsequently, they were predicted by the mtk-QSBER model. Regarding the reliability of the predictions for the designed molecules, we applied two criteria. First, our decision of classifying each new molecule as positive was based on a consensus prediction. In this sense, each designed molecule was predicted three times against each of the 230 combinations of the elements me, bt, ai, and tm (already reported for the calculation of the fragment contributions). The three times correspond to the three degrees of reliability of the experimental information, i.e., autocuration (pc = 0.55), intermediate (pc = 0.75), and expert (pc = 1). If in at least 2 out of 3 possible times a designed molecule was predicted as positive against a defined combination of me, bt, ai, and tm, then the molecule was annotated as positive. Otherwise, the molecule was assigned as negative. At the end, each of the six new molecules was predicted in 230 × 3 = 690 experimental conditions. That yielded a total of 690 × 6 molecules = 4140 predictions involving both anti-HCV activity and ADMET profiles. Second, each prediction performed for any of the designed molecules was considered reliable if it was within the applicability domain of the mtk-QSBER model. Here, the applicability domain was defined according to the bounding box approach,57 that is, a multi-dimensional hyper-rectangle defined on the basis of the maximum and minimum values of each molecular descriptor used to build the mtk-QSBER model (Table 5). Notice that the applicability domain was determined only with the compounds present in the training set. The results of the virtual predictions can be found in the Supplementary Information 5. The six designed molecules were predicted as positive against most of the 690 experimental conditions. In general, the screening of these new molecules suggests that they should have potent inhibitory activity 10

ACS Paragon Plus Environment

Page 10 of 23

Page 11 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

against all the eight HCV strains present in our dataset, as well as very good permeability, and a relatively low protein plasma binding. From the predictions, one can also assume that the molecules may be weak inhibitors almost all the metabolizing proteins (30 in total) reported in this study, namely, cytochromes, transferases, and transporters. The only exceptions were registered for cytochrome P450 2A6, where the molecules AHV-3 and AHV-6 were predicted as negative in 9 out of 9 times against the aforementioned enzyme (Supplementary Information 5). Therefore, these two molecules may exhibit IC50 (nM)cyp < 5000 nM according to the cutoff depicted in Table 1. At the same time, the molecule AHV-1 was predicted as positive in 7 out of 9 times, while for AHV-5, the same result was achieved in 5 out of 9 times. Thus, both molecules can still be viewed as potentially weak cytochrome P450 2A6 inhibitors. Molecules AHV-2 and AHV-4 were classified as positive in 9 out of 9 times, which suggests that they will have a better chance of being confirmed as weak inhibitors than the other designed molecules. Table 5. Applicability domain for the mtk-QSBER (bounding box approach). Descriptor DBq3(AR)m e DBq1(AW)bt DBq7(HYD)ai DBq1(PSA)ai DBq6(AW)tm

Minimum

Maximum

-3021.68 15168.29 -91.01 337.53 -21870.7 139178.21 -1171.65 8873.19 -31933.27 206261.93

Other prediction results indicate that all the new molecules will exhibit low clearance rate and large half lives. From a toxicological point of view, at the molecular level, the mtk-QSBER model predicted the six designed molecules to be weak inhibitors of the 11 ion channels (including HERG) reported in this study. In convergence with this result, at the cellular level, the molecules were predicted to exhibit very low toxic effects against all the mammalian cells (14 in total). Altogether, the toxicity predictions recommend that the six molecules designed here may lack various serious side effects such as cardiotoxicity, hepatotoxicity, and nephrotoxicity, among others. On the other hand, and also in terms of reliability, only 20 out of 4140 predictions mentioned above fell outside the applicability domain of the mtk-QSBER model. Such predictions were obtained for the targets labeled as “sulfonylurea receptors; K-ATP channels”, “sodium channel alpha subunits; brain (types I, II, III)”, and “dopamine transporter” only against the molecules AHV-1, AHV-3, AHV-4, and AHV-5. For each of them, the same following result was obtained: reliable predictions were reported in 4 out of 6 times against each of the targets “sulfonylurea receptors; K-ATP channels” and “sodium channel alpha subunits; brain (types I, II, III)”, while for the target “dopamine transporter”, in 23 out of 24 times the predictions were considered reliable. Finally, we would like to point out that the six designed molecules were not reported in our dataset. In fact, we also searched for similar molecules in drugs/chemicals databases such as CHEMBL,30 PubChem,58 ZINC,59 and eMolecules.60 In all of them, we applied a cutoff ≥ 80% as the criterion in the search for molecular similarity. No results were found. Assessing the druglikeness of the designed molecules One of the most important stages in drug discovery is the assessment of the druglikeness of a molecule, which is strongly related to its oral bioavailability. In many drug discovery projects in both academia and pharmaceutical industry, a traditional method to evaluate druglikeness is to check compliance with the Lipinski’s rule of five.61 According to this rule, and from the perspective of a medicinal chemist, any molecule considered to be orally administered should possess the following features: 1) no more than 5 hydrogen bond donors (the total number of nitrogen–hydrogen and oxygen–hydrogen bonds); 2) no more than 10 hydrogen bond acceptors (all nitrogen or oxygen atoms); 3) a molecular mass less than 500 daltons; 4) a logarithm of the octanol-water partition coefficient (logP) not greater than 5 or Moriguchi’s logP not greater than 4.15. A second variant62 of the aforementioned rule is logP in −0.4 to +5.6 range, a molar refractivity from 40 to 130, a molecular weight from 180 to 500, a number of atoms from 20 to 70, and a polar surface area no greater than 140 Ǻ2. A third variant63 has criticized Lipinski’s rule of five regarding the cutoff of the molecular mass (500 daltons), suggesting that only molecules with 10 or fewer rotatable bonds and polar surface area (only nitrogen and oxygen atoms considered) equal to or less than 140 Å2 (or 12 or fewer H-bond donors and acceptors) are predicted to have good oral bioavailability. In order to check the druglikeness of the designed molecules, we calculated all the aforementioned 11

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 23

properties using the software DRAGON.64 All the properties appear depicted in Table 6. In the specific case of the property logP, we used two different approaches, which are available in DRAGON (MlogP Moriguchi’s logP),64 and ChemDraw Ultra (ClogP).65 According to Table 6, there are no violations of the Lipinski’s rule of five or any of its variants. Even in the specific case of the property logP calculated by different approaches, one can see that despite the divergences in the values, all are in the interval 1.624.08, which is in accordance with Lipinski’s rule of five. Table 6. Different molecular properties of the designed molecules. IDa,b

MW

nAT

RBN

nHDon

nHAcc

AMR

TPSA(Tot)

MlogP

ClogP

AHV-1

414.96

47

8

3

6

105.39

116.13

2.61

2.88

AHV-2

431.91

44

7

3

7

103.52

132.03

2.92

2.97

AHV-3

380.87

39

8

1

8

87.11

113.95

2.95

2.24

AHV-4

399.43

42

7

2

8

94.78

125.72

1.88

3.32

AHV-5

415.91

43

7

2

6

100.09

103.96

4.08

2.85

AHV-6 430.85 44 10 0 7 92.06 112.5 1.62 1.72 The symbols refer to the following properties: MW – molecular weight, nAT – number of atoms, RBN – number of rotatable bonds, nHDon – number of hydrogen bond donors, nHAcc – number of hydrogen bond acceptors, AMR – molar refractivity, TPSA(Tot) – polar surface area based on fragments containing nitrogen, oxygen, phosphorus, and sulfur, MlogP – logarithm of the partition coefficient (octanol/water) according to the Moriguchi’s approach available in the software DRAGON, ClogP – logarithm of the partition coefficient (octanol/water) calculated by the program ChemDraw Ultra v8.0. b For the calculation of the property nHAcc, fluorine atoms are also considered as hydrogen bond acceptors. a

Conclusions The search for new anti-hepatitis drugs remains a challenging scientific area. Due to the fast growing body of high-quality available experimental data, computational models should be employed beyond the classical objective of screening databases. When possible, more mechanistic and phenomenological information should be extracted from the different molecular descriptors in order to create more biologically meaningful models. Special attention should be paid to the fact that in silico models may be used as tools to perform true rational design of new therapeutic agents. In this sense, the mtk-QSBER model derived here can be viewed as an encouraging alternative in antiviral research for the discovery of potentially promising anti-hepatitis C leads. The very good statistical quality and predictive power of the mtk-QSBER model, along with the physicochemical interpretations of the molecular descriptors, and the use of the fragments contributions allowed us to design six new molecules that were predicted to display potent anti-HCV and suitable in vitro ADMET properties. The six new molecules may be used for the generation of focused chemical libraries. The present mtk-QSBER model creates the foundation to speed up the virtual design and screening of versatile and highly efficacious anti-hepatitis C drugs, serving as a complimentary tool of powerful techniques such as HTS in antiviral research. Acknowledgments This work received financial support from Fundação para a Ciência e a Tecnologia (FCT/MEC) through national funds, and co-financed by the European Union (FEDER funds) under the Partnership Agreement PT2020, through projects UID/QUI/50006/2013, POCI/01/0145/FEDER/007265, and NORTE-01-0145FEDER-000011 (LAQV@REQUIMTE). To all financing sources, the authors are greatly indebted. Supporting Information Chemical and biological data, and averages (Supplementary Information 1.xlsx). Input file and classification results (Supplementary Information 2.xlsx). %CCC(cj) values (Supplementary Information 3.xlsx). Fragment contributions (Supplementary Information 4.xlsx). Virtual screening performed for the designed molecules (Supplementary Information 5.xlsx).

12

ACS Paragon Plus Environment

Page 13 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

References (1) Wree, A.; Broderick, L.; Canbay, A.; Hoffman, H. M.; Feldstein, A. E. From NAFLD to NASH to cirrhosis-new insights into disease mechanisms. Nat. Rev. Gastroenterol. Hepatol. 2013, 10, 627-636. (2) Bernal, W.; Auzinger, G.; Dhawan, A.; Wendon, J. Acute liver failure. Lancet 2010, 376, 190-201. (3) Forner, A.; Llovet, J. M.; Bruix, J. Hepatocellular carcinoma. Lancet 2012, 379, 1245-1255. (4) Mohd Hanafiah, K.; Groeger, J.; Flaxman, A. D.; Wiersma, S. T. Global epidemiology of hepatitis C virus infection: new estimates of age-specific antibody to HCV seroprevalence. Hepatology 2013, 57, 1333-1342. (5) Mesalam, A. A.; Vercauteren, K.; Meuleman, P. Mouse systems to model hepatitis C virus treatment and associated resistance. Viruses 2016, 8, 176. (6) Negro, F. Adverse effects of drugs in the treatment of viral hepatitis. Best Pract. Res. Clin. Gastroenterol. 2010, 24, 183-192. (7) Paul, S. M.; Mytelka, D. S.; Dunwiddie, C. T.; Persinger, C. C.; Munos, B. H.; Lindborg, S. R.; Schacht, A. L. How to improve R&D productivity: the pharmaceutical industry's grand challenge. Nat. Rev. Drug Discov. 2010, 9, 203-214. (8) Jahnke, W.; Erlanson, D. A. Fragment-based approaches in drug discovery; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2006. (9) Liao, C.; Sitzmann, M.; Pugliese, A.; Nicklaus, M. C. Software and resources for computational medicinal chemistry. Future Med. Chem. 2011, 3, 1057-1085. (10) Vrontaki, E.; Melagraki, G.; Mavromoustakos, T.; Afantitis, A. Searching for anthranilic acid-based thumb pocket 2 HCV NS5B polymerase inhibitors through a combination of molecular docking, 3DQSAR and virtual screening. J. Enzyme Inhib. Med. Chem. 2016, 31, 38-52. (11) Yu, H.; Fang, Y.; Lu, X.; Liu, Y.; Zhang, H. Combined 3D-QSAR, molecular docking, molecular dynamics simulation, and binding free energy calculation studies on the 5-hydroxy-2H-pyridazin-3-one derivatives as HCV NS5B polymerase inhibitors. Chem. Biol. Drug Des. 2014, 83, 89-105. (12) Worachartcheewan, A.; Prachayasittikul, V.; Toropova, A. P.; Toropov, A. A.; Nantasenamat, C. Large-scale structure-activity relationship study of hepatitis C virus NS5B polymerase inhibition using SMILES-based descriptors. Mol. Divers. 2015, 19, 955-964. (13) Speck-Planche, A.; Cordeiro, M. N. D. S. Computer-aided drug design methodologies toward the design of anti-hepatitis C agents. Curr. Top. Med. Chem. 2012, 12, 802-813. (14) Therese, P. J.; Manvar, D.; Kondepudi, S.; Battu, M. B.; Sriram, D.; Basu, A.; Yogeeswari, P.; Kaushik-Basu, N. Multiple e-pharmacophore modeling, 3D-QSAR, and high-throughput virtual screening of hepatitis C virus NS5B polymerase inhibitors. J. Chem. Inf. Model. 2014, 54, 539-552. (15) Wei, Y.; Li, J.; Qing, J.; Huang, M.; Wu, M.; Gao, F.; Li, D.; Hong, Z.; Kong, L.; Huang, W.; Lin, J. Discovery of Novel Hepatitis C Virus NS5B Polymerase Inhibitors by Combining Random Forest, Multiple e-Pharmacophore Modeling and Docking. PLoS One 2016, 11, e0148181. (16) Qureshi, A.; Kaur, G.; Kumar, M. AVCpred: an integrated web server for prediction and design of antiviral compounds. Chem. Biol. Drug Des. 2017, 89, 74-83. (17) Romero-Duran, F. J.; Alonso, N.; Yanez, M.; Caamano, O.; Garcia-Mera, X.; Gonzalez-Diaz, H. Brain-inspired cheminformatics of drug-target brain interactome, synthesis, and assay of TVP1022 derivatives. Neuropharmacology 2016, 103, 270-278. (18) Tenorio-Borroto, E.; Penuelas-Rivas, C. G.; Vasquez-Chagoyan, J. C.; Castanedo, N.; Prado-Prado, F. J.; Garcia-Mera, X.; Gonzalez-Diaz, H. Model for high-throughput screening of drug immunotoxicity Study of the anti-microbial G1 over peritoneal macrophages using flow cytometry. Eur. J. Med. Chem. 2014, 72, 206-220. (19) Speck-Planche, A.; Kleandrova, V. V.; Ruso, J. M.; Cordeiro, M. N. D. S. First multitarget chemobioinformatic model to enable the discovery of antibacterial peptides against multiple Gram-positive pathogens. J. Chem. Inf. Model. 2016, 56, 588-598. (20) Speck-Planche, A.; Cordeiro, M. N. D. S. Enabling virtual screening of potent and safer antimicrobial agents against noma: mtk-QSBER model for simultaneous prediction of antibacterial activities and ADMET properties. Mini Rev. Med. Chem. 2015, 15, 194-202. (21) Speck-Planche, A.; Cordeiro, M. N. D. S. Chemoinformatics for medicinal chemistry: in silico model to enable the discovery of potent and safer anti-cocci agents. Future Med. Chem. 2014, 6, 20132028. (22) Zhang, W.; Pei, J.; Lai, L. Computational Multitarget Drug Design. J. Chem. Inf. Model. 2017, 57, 403-412. (23) Hodos, R. A.; Kidd, B. A.; Shameer, K.; Readhead, B. P.; Dudley, J. T. In silico methods for drug repurposing and pharmacology. Wiley Interdiscip. Rev. Syst. Biol. Med. 2016, 8, 186-210. 13

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 23

(24) Speck-Planche, A.; Cordeiro, M. N. D. S. Multitasking models for quantitative structure-biological effect relationships: current status and future perspectives to speed up drug discovery. Expert Opin. Drug Discov. 2015, 10, 245-256. (25) Herrera-Ibata, D. M.; Pazos, A.; Orbegozo-Medina, R. A.; Romero-Duran, F. J.; Gonzalez-Diaz, H. Mapping chemical structure-activity information of HAART-drug cocktails over complex networks of AIDS epidemiology and socioeconomic data of U.S. counties. Biosystems 2015, 132-133, 20-34. (26) Herrera-Ibata, D. M.; Pazos, A.; Orbegozo-Medina, R. A.; Gonzalez-Diaz, H. Mapping networks of anti-HIV drug cocktails vs. AIDS epidemiology in the US counties. Chemometr. Intell. Lab. Syst. 2014, 138, 161-170. (27) Gonzalez-Diaz, H.; Herrera-Ibata, D. M.; Duardo-Sanchez, A.; Munteanu, C. R.; Orbegozo-Medina, R. A.; Pazos, A. ANN multiscale model of anti-HIV drugs activity vs AIDS prevalence in the US at county level based on information indices of molecular graphs and social networks. J. Chem. Inf. Model. 2014, 54, 744-755. (28) Kleandrova, V. V.; Speck Planche, A. Multitasking model for computer-aided design and virtual screening of compounds with high anti-HIV activity and desirable ADMET properties. In Multi-Scale Approaches in Drug Discovery: From Empirical Knowledge to In Silico Experiments and Back, 1st ed.; Speck-Planche, A., Ed. Elsevier: Oxford, UK, 2017; pp 55-81. (29) Speck-Planche, A.; Cordeiro, M. N. D. S. Simultaneous modeling of antimycobacterial activities and ADMET profiles: a chemoinformatic approach to medicinal chemistry. Curr. Top. Med. Chem. 2013, 13, 1656-1665. (30) Gaulton, A.; Bellis, L. J.; Bento, A. P.; Chambers, J.; Davies, M.; Hersey, A.; Light, Y.; McGlinchey, S.; Michalovich, D.; Al-Lazikani, B.; Overington, J. P. ChEMBL: a large-scale bioactivity database for drug discovery. Nucleic Acids Res. 2012, 40, D1100-1107. (31) Mok, N. Y.; Brenk, R. Mining the ChEMBL database: an efficient chemoinformatics workflow for assembling an ion channel-focused screening library. J. Chem. Inf. Model. 2011, 51, 2449-2454. (32) Overington, J. ChEMBL. An interview with John Overington, team leader, chemogenomics at the European Bioinformatics Institute Outstation of the European Molecular Biology Laboratory (EMBLEBI). Interview by Wendy A. Warr. J. Comput. Aided Mol. Des. 2009, 23, 195-198. (33) Tenorio-Borroto, E.; Garcia-Mera, X.; Penuelas-Rivas, C. G.; Vasquez-Chagoyan, J. C.; PradoPrado, F. J.; Castanedo, N.; Gonzalez-Diaz, H. Entropy model for multiplex drug-target interaction endpoints of drug immunotoxicity. Curr. Top. Med. Chem. 2013, 13, 1636-1649. (34) ChemAxon Standardizer (Tool for structure canonicalization and transformation), JChem, v15.11.16.0; Budapest, Hungary, 1998-2016. (35) Valdés-Martini, J. R.; García-Jacas, C. R.; Marrero-Ponce, Y.; Silveira Vaz ‘d Almeida, Y.; Morell, C. QUBILs-MAS: Free software for molecular descriptors calculator from quadratic, bilinear and linear maps based on graph-theoretic electronic-density matrices and atomic weightings, v1.0; CAMD-BIR Unit, CENDA registration number: 2373-2012: Villa Clara, 2012. (36) Marrero-Ponce, Y.; Torrens, F.; Alvarado, Y. J.; Rotondo, R. Bond-based global and local (bond, group and bond-type) quadratic indices and their applications to computer-aided molecular design. 1. QSPR studies of diverse sets of organic chemicals. J. Comput. Aided Mol. Des. 2006, 20, 685-701. (37) Meneses-Marcel, A.; Rivera-Borroto, O. M.; Marrero-Ponce, Y.; Montero, A.; Machado Tugores, Y.; Escario, J. A.; Gomez Barrio, A.; Montero Pereira, D.; Nogal, J. J.; Kouznetsov, V. V.; Ochoa Puentes, C.; Bohorquez, A. R.; Grau, R.; Torrens, F.; Ibarra-Velarde, F.; Aran, V. J. New antitrichomonal drug-like chemicals selected by bond (edge)-based TOMOCOMD-CARDD descriptors. J. Biomol. Screen. 2008, 13, 785-794. (38) Casanola-Martin, G. M.; Marrero-Ponce, Y.; Khan, M. T.; Khan, S. B.; Torrens, F.; Perez-Jimenez, F.; Rescigno, A.; Abad, C. Bond-based 2D quadratic fingerprints in QSAR studies: virtual and in vitro tyrosinase inhibitory activity elucidation. Chem. Biol. Drug Des. 2010, 76, 538-545. (39) Castillo-Garit, J. A.; Vega, M. C.; Rolon, M.; Marrero-Ponce, Y.; Kouznetsov, V. V.; Torres, D. F.; Gomez-Barrio, A.; Bello, A. A.; Montero, A.; Torrens, F.; Perez-Gimenez, F. Computational discovery of novel trypanosomicidal drug-like chemicals by using bond-based non-stochastic and stochastic quadratic maps and linear discriminant analysis. Eur. J. Pharm. Sci. 2010, 39, 30-36. (40) Kleandrova, V. V.; Ruso, J. M.; Speck-Planche, A.; Dias Soeiro Cordeiro, M. N. Enabling the Discovery and Virtual Screening of Potent and Safe Antimicrobial Peptides. Simultaneous Prediction of Antibacterial Activity and Cytotoxicity. ACS Comb. Sci. 2016, 18, 490-498. (41) Hill, T.; Lewicki, P. STATISTICS Methods and Applications. A Comprehensive Reference for Science, Industry and Data Mining; StatSoft: Tulsa, 2006. (42) Statsoft-Team STATISTICA. Data analysis software system, v6.0; Tulsa, 2001. 14

ACS Paragon Plus Environment

Page 15 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

(43) Cronin, M. T. D.; Schultz, T. W. Pitfalls in QSAR. Theochem-J. Mol. Struc. 2003, 622, 39-51. (44) Pearson, K. Notes on regression and inheritance in the case of two parents. Proc. R. Soc. Lond. 1895, 58, 240-242. (45) Chauhan, J. S.; Dhanda, S. K.; Singla, D.; Agarwal, S. M.; Raghava, G. P. QSAR-based models for designing quinazoline/imidazothiazoles/pyrazolopyrimidines based inhibitors against wild and mutant EGFR. PLoS One 2014, 9, e101079. (46) Pearson, K. On lines and planes of closest fit to systems of points in space. Philos. Mag. 1901, 2, 559-572. (47) Hotelling, H. Analysis of a complex of statistical variables into principal components. J. Educ. Psychol. 1933, 24, 417-441. (48) Hotelling, H. Relations between two sets of variates. Biometrika 1936, 28, 321-377. (49) Randic, M. Orthogonal molecular descriptors. New J. Chem. 1991, 15, 517-525. (50) Todeschini, R.; Consonni, V. Handbook of Molecular Descriptors; WILEY-VCH Verlag GmbH: Weinheim, New York, Chichester, Brisbane, Singapore, Toronto, 2000. (51) Baskin, I. I.; Skvortsova, M. I.; Stankevich, I. V.; Zefirov, N. S. On the basis of invariants of labeled molecular graphs. J. Chem. Inf. Comput. Sci. 1995, 35, 527-531. (52) Baskin, I.; Varnek, A. Fragment descriptors in SAR/QSAR/QSPR studies, molecular similarity analysis and in virtual screening. In Chemoinformatics approaches to virtual screening, Varnek, A.; Tropsha, A., Eds; Royal Society of Chemistry: Cambridge, 2008; pp 1-43. (53) Estrada, E.; Uriarte, E.; Montero, A.; Teijeira, M.; Santana, L.; De Clercq, E. A novel approach for the virtual screening and rational design of anticancer compounds. J. Med. Chem. 2000, 43, 1975-1985. (54) Toplak, M.; Mocnik, R.; Polajnar, M.; Bosnic, Z.; Carlsson, L.; Hasselgren, C.; Demsar, J.; Boyer, S.; Zupan, B.; Stalring, J. Assessment of machine learning reliability methods for quantifying the applicability domain of QSAR regression models. J. Chem. Inf. Model. 2014, 54, 431-441. (55) Carrio, P.; Pinto, M.; Ecker, G.; Sanz, F.; Pastor, M. Applicability Domain ANalysis (ADAN): a robust method for assessing the reliability of drug property predictions. J. Chem. Inf. Model. 2014, 54, 1500-1511. (56) Gaspar, H. A.; Marcou, G.; Horvath, D.; Arault, A.; Lozano, S.; Vayer, P.; Varnek, A. Generative topographic mapping-based classification models and their applicability domain: application to the biopharmaceutics Drug Disposition Classification System (BDDCS). J. Chem. Inf. Model. 2013, 53, 3318-3325. (57) Sahigara, F.; Mansouri, K.; Ballabio, D.; Mauri, A.; Consonni, V.; Todeschini, R. Comparison of different approaches to define the applicability domain of QSAR models. Molecules 2012, 17, 4791-4810. (58) Kaiser, J. Science resources. Chemists want NIH to curtail database. Science 2005, 308, 774. (59) Irwin, J. J.; Shoichet, B. K. ZINC--a free database of commercially available compounds for virtual screening. J. Chem. Inf. Model. 2005, 45, 177-182. (60) Gubernator, K.; James, C. A.; Gubernator, N. eMolecules. https://www.emolecules.com/ (February 10, 2017). (61) Lipinski, C. A.; Lombardo, F.; Dominy, B. W.; Feeney, P. J. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv. Drug Deliv. Rev. 2001, 46, 3-26. (62) Ghose, A. K.; Viswanadhan, V. N.; Wendoloski, J. J. A knowledge-based approach in designing combinatorial or medicinal chemistry libraries for drug discovery. 1. A qualitative and quantitative characterization of known drug databases. J. Comb. Chem. 1999, 1, 55-68. (63) Veber, D. F.; Johnson, S. R.; Cheng, H. Y.; Smith, B. R.; Ward, K. W.; Kopple, K. D. Molecular properties that influence the oral bioavailability of drug candidates. J. Med. Chem. 2002, 45, 2615-2623. (64) Todeschini, R.; Consonni, V.; Mauri, A.; Pavan, M. DRAGON for Windows (Software for Molecular Descriptor Calculations), v5.3; Milano Chemometrics and QSAR Research Group: Milano, 2005. (65) CambridgeSoft ChemDraw Ultra, v8.0; Cambridge, MA, 2003.

15

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 23

Captions for figures Figure 1. Diverse stages involved in the creation of the mtk-QSBER model. Figure 2. AUROC curves. Figure 3. Standardized coefficients used as measures of the relative influences of the molecular descriptors. Figure 4. Fragments selected for the calculation of the quantitative contributions. Figure 5. Structures of the new molecules designed by using the mtk-QSBER model.

16

ACS Paragon Plus Environment

Page 17 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

Graphical Abstract (tif format) 28x8mm (600 x 600 DPI)

ACS Paragon Plus Environment

1 2CHEMBL DATABASE 3 4 5 6 7

SMILES

Ki, IC , Combinatorial Science ACS 50 EC50, PPB, CL, t1/2, CC50

ACS Paragon Plus Environment Proteins

Hepatitis C Virus Mammalian cells

Virtual design

Page 18 ofand 23 of potent

safe anti-hepatitis C agents

Page 19 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

Figure 1. Diverse stages involved in the creation of the mtk-QSBER model. 21x18mm (600 x 600 DPI)

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. AUROC curves. 14x11mm (600 x 600 DPI)

ACS Paragon Plus Environment

Page 20 of 23

Page 21 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

Figure 3. Standardized coefficients used as measures of the relative influences of the molecular descriptors. 11x6mm (600 x 600 DPI)

ACS Paragon Plus Environment

ACS Combinatorial Science

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Fragments selected for the calculation of the quantitative contributions. 203x226mm (600 x 600 DPI)

ACS Paragon Plus Environment

Page 22 of 23

Page 23 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Combinatorial Science

Figure 5. Structures of the new molecules designed by using the mtk-QSBER model. 194x214mm (600 x 600 DPI)

ACS Paragon Plus Environment