Antiprotozoal Nitazoxanide Derivatives: Synthesis, Bioassays and QSAR Study Combined with Docking for Mechanistic Insight.

In view of the serious health problems concerning infectious diseases in heavily populated areas, we followed the strategy of lead compound diversification to evaluate the near-by chemical space for new organic compounds. To this end, twenty derivatives of nitazoxanide (NTZ) were synthesized and tested for activity against Entamoeba histolytica parasites. To ensure drug-likeliness and activity relatedness of the new compounds, the synthetic work was assisted by a quantitative structure-activity relationships study (QSAR). Many of the inherent downsides - well-known to QSAR practitioners - we circumvented thanks to workarounds which we proposed in prior QSAR publication. To gain further mechanistic insight on a molecular level, ligand-enzyme docking simulations were carried out since NTZ is known to inhibit the protozoal pyruvate ferredoxin oxidoreductase (PFOR) enzyme as its biomolecular target.

INTRODUCTION

In the 1930s, chemists were beginning to explore the effects of structural modifications on the kinetic processes of chemical reactions, resulting in the birth of physical organic chemistry [1]. Decades before, other scientists observed structure-activity relationships, e.g. ethers or alcohols (Cros, “Action de l'alcohol amylique sur l'organisme”, University of Strasbourg, France, 1863) showed a correlation between lipo-solubility and toxicity. Around 1900, Meyer and Overton, independently, established the linear dependency between the narcotic action and water / oil partitioning (Meyer, 1899; Overton, 1901 in [2]). One of the seminal works at an early stage of Quantitative Structure-Activity Relationships (QSAR) was the study on electronic effects of benzoic acids substituents and the descriptors were named after its inventor Hammet σ constants (Hammet, 1937; Hammet, 1940 in [2]). After the 1960s, QSAR modeling was carried out using multiple linear regression (MLR) on numeric independent variables describing structural features (descriptors) by Hansch and Fujita (Hansch, Fujita, 1964 in [2]) under the assumption that an “inherent association between chemical structure and biological activity” exists [3]. Not more than twenty physicochemical descriptors are commonly used in QSAR studies descriptors [4], like octanol/water partition coefficient (log P), the Hammet σ constant, and the fragmental lipophilicity parameter π. Others were derived from quantum chemical calculations, namely: dipole moment, partial charges, HOMO/LUMO energies (Kubinyi, 1993 in [2]) or were based on molecular graph theory and topology concepts by Wiener (Wiener, 1947 in [2]). Kier and Hall (Kier, Hall, Murray, 1975 in [2]) and Randic (Randic, 1975; Kier, Hall, 1976; Kier, Hall, 1986 in [2]).

QSAR studies were applied on many occasions and recently reviewed [5-7]. QSAR studies also predicted antiprotozoal potency of pentamidines. The authors generated their own descriptors and used them in the final multiple linear regression equations. The comparison to our descriptors is delayed, awaiting the free Desmol software [8]. Drug screens from natural and synthetic sources have focused on antiprotozoal activities as part of the Drugs for Neglected Diseases Initiative [9-13]. A kind of knowledge-based approach is the combination of activity data base and QSAR prediction, called multi-species complex networks for antiparasite drugs [14]. As early as the 1980s nitro-aromatic antiparasitic drugs were already well-known, e.g. metronidazole [15]. Structurally related compounds to our nitrothiazole scaffold are also reported, e.g. benznidazole, nifurtimox and megazol 5-nitroimidazole derivative highly active in vitro and in vivo against parasites like Trypanosoma cruzi, or thiadiazines [16, 17].

With growing population worldwide, infectious diseases are considered a serious health problem. In particular, amoebiasis is a protozoal infection caused by a living organism called Entamoeba histolytica. Upon contaminated food uptake or other oral expositions, the clinical symptoms of patients include severe dysenteries, mucosal ulcers and even peritonitis, amoebic granulomas and/or fulminant colitis [18-20].

In 1974, Rossignol and Cavier conducted synthesis of nitazoxanide (NTZ) which was patented two years later [21]. Nowadays the drug is employed for a plethora of treatments, such as: anaerobic bacteria (Helicobacter pylori, Clostridium difficile); protozoa (, Balantidium coli, Giardia intestinalis, Trichomonas vaginalis); helminthes (E. vermicularis, A. lumbricoides, Strongyloides stercoralis, T. trichuria, Taenia spp, H.nana, Fasciola hepática) in addition to certain viruses (rotavirus, hepatitis B and C and influenza) [22-26].

The common sense indicates that the biological activity of NTZ relies on the presence of the 5-nitrothiazole ring as the core substructure, i.e. the essential scaffold (pharmacophore) which is kept constant during chemical derivatization efforts (Figs. , ) [27]. Deacetylation of NTZ produces tizoxanide (TIZ). It is a known active metabolite and included in the study (Tables and ) [27-30].

In 2011, Navarrete-Vazquez et al. synthesized and characterized two benzologues of nitazoxanide and tizoxanide [27], with proven antiprotozoal activity. After that, new compounds were designed and synthesized (Fig. ), and clustered into two series containing 11 thiazole analogues and 9 benzothiazoles as benzologues. They all (n = 20) were tested for their antiprotozoal activities (IC50) against Entamoeba histolytica [personal communication, 2011]. The molecular structures are displayed in Figs. ), Tables and . Moreover, Hoffman suggested NTZ acts as an inhibitor of the pyruvate ferredoxin oxidoreductase (PFOR) enzyme as a possible mechanism of action [31, 32, personal communication 2012].

Pharmacological reports relate the antiprotozoal activity of NTZ to an anaerobic microorganism-specific enzyme as its drug target, namely PFOR [31, 32, 35-37]. The physiological role of PFOR is known as an oxidative decarboxylation of pyruvate to form acetyl-CoA and carbon dioxide, paralleled by the reduction of ferredoxin or flavodoxin in the presence of the coenzyme Thiamine Pyrophosphate (TPP, Vitamin B1) [32, 37-39]. One possible explanation for the essential step in the anaerobic energy metabolism was proposed in that the antiprotozoal drug (candidates) inhibits the electron transfer reaction in a PFOR-dependent fashion [24, 30, 32]. Another possible pathway was suggested to explain the antiprotozoal activity [40]. P. S. Hoffman et al. proposed that NTZ directly interacts with TPP, what is a unique case of drug targeting a coenzyme. They observed NTZ was not replacing TPP from the active site at all in their Helicobacter pylori and Campylobacter jejuni assays [32].

On occasion of practicing a full-size antiprotozoal drug development cycle including in silico design, syntheses and in vitro bioassays, we also inspected the quantitative structure-activity relationships (QSAR) models and analyzed their equations in order to identify the pitfalls encountered therein. In an additional step we tried to circumvent them if possible according to our 2009 report on a collection of frequent problems and difficulties met in QSAR studies [3]. Intriguingly, A. Hopfinger et al. successfully applied predictive 3D- and 4D-QSAR techniques to two series of pyridinyl- and pyrimidinyl-imidazoles [41, 42], so we decided to apply higher dimensional QSAR protocols for our NTZ derivatives, too.

For more details concerning modeling, bioassays and mathematical model, including leave one out cross-validation procedures, see the Supplementary Material.

MATERIALS AND METHODS

Organic Synthesis and Characterization of New Compounds

Triethylamine (1.2 equiv.) was added to a solution of 2-amino-5-nitrothiazole or 2-amino-6-nitrobenzo[d]thiazole (0.0015 mol) in dicloromethane. Fifteen minutes of stirring was necessary at 5°C. To this mixture a solution of acylchloride (adequately substituted) in dichloromethane was added. Under continuous stirring at room temperature during several hours the reaction process was monitored by thin layer chromatography between 4 to 24 hours until the transformation was complete. After distillation of dichloromethane was achieved, the resulting residues were neutralized with a saturated NaHCO3 solution. Prior to washing with brine (3×20 mL) the aqueous layer was extracted with ethyl acetate (3×15 mL). Ethyl acetate was removed in vacuum. Then the precipitated solids were either recrystallized from a mixture of solvents or purified by column chromatography.

QSAR Modeling and Molecular Docking

Upon three-dimensional (3D) model generations [43-46], their molecular geometries adopted a conformation taken from the crystal structure of NTZ [47, 48]. Prior to retrieving the crystal structure of NTZ, we correctly predicted (15°) the out of coplanarity tilt (18°) between the 5- and 6-membered rings of NTZ.

Prior to docking of NTZ into the protozoal enzyme PFOR, the missing structure of the target enzyme PFOR from Entamoeba histolytica (RCSB Protein Data Bank) [49] was built from a similar PFOR crystal structure of Desulfovibrio africanus as three-dimensional template (PDB code: 2C42 [50]) [51, 52]. Automated flexible ligand docking was applied to the PFOR model of E. histolytica as rigid target receptor under Auto-Dock 4.2 & AutoDockTools [53, 54]. Adjacent to the binding cleft we kept the prosthetic group ferrodoxin and a magnesium cation to give the site its natural form.

Parts of our multicenter molecular modeling contributions were carried out on a Linux workstation using Schrodinger/CANVAS software [55]. To obtain QSAR equations by stepwise forward MLR method, an in-house program was developed using C language on a Windows platform.

Various 2D topological, structural, constitutional and physicochemical descriptors were calculated using CANVAS yielding a set of 757 descriptors [55]. In addition, a binary variable which has values 0 or 1 for indicating absence or presence of benzothiazole scaffold respectively was included in the study. The data preprocessing of independent variables was performed by removing the descriptors common to more than 90% of the compounds and with a homogeneity index < 20 (i.e. deviating more than 80% from the ideal homogenous distribution). The homogeneity index quantitatively indicates the degree of uniformity of a variable (descriptor) distribution within its range [56]. Data preprocessing led to a significant reduction of descriptors viz. 217 descriptors for further study. Furthermore for building multiple linear regression based QSAR models auto-scaled descriptors were utilized.

For the statistical procedure we followed the reasoning about internal and external validation processes by M.N. Noolvi et al. and applied their equations appropriately, besides alternative solutions [57, 58]. On the other side, the drawbacks and shortcomings concerning cross-validation procedures were published in a seminal work by D. Krstajic and coworkers [59]. I. Tetko (see below for Ochem online QSAR server [60]) recommended that the final step should always be the validation of hitherto unseen compounds and done only once (private communications, Helmholtz-Zentrum Munich, Germany, 2013, 2014). We followed this tenet wherever the software allowed us to maintain all 22 molecules together for QSAR modeling input, while Raptor was the sole tool which requires an a prior manual division into a test and training set [61, 62]. Hence our data set was only divided into two groups for Raptor application. To this end, the Raptor training set contained 18 compounds to generate the 5D-QSAR model complemented by a test set with 4 hand-selected compounds for the external (final) model validation. The molecules labeled B10, B30, T04 and T18 were manually selected because they embrace almost all the chemical features of the training set. The underlying trade-off meant to reflect a minimum number of molecules in the test set, while maximizing the chemical representativeness in that test set. The leave one out (LOO, Q2) method was carried out for internal/cross-validation of the models. For validation of either internal/cross or external models, values of either Q2 or pred R2, and other statistical parameters (Q2F1, Q2F2, Q2F3, r2m) were calculated using formulae reported in the literature and references cited therein [57, 63-70].

While Q2 values describe the internal stability, the pred R2 values reflects the predictive power of a model [57]. Moreover, to the influence of chance correlations on the model building of the training set was evaluated by Y-randomization (randomization of the biological activity) test and ‘Z score’ as reported in literature [3, 71]. The Y-randomization test was performed to ascertain that the obtained statistical significance of the QSAR model was not due to chance correlation or overfitting.

With another tool box, we created additional QSAR linear equations using widely accepted concepts [44, 72-75]. The acidic dissociation constants (pKa) were estimated in the first place based on the known functional fragments and later determined by two online resources (Marvin, SPARC) [76, 77]. A higher-dimensional QSAR approach included a 4D-QSAR [45, 46] as well as 5D-QSAR study [61, 62]. For external validation of the QSAR equations, newly, all molecules were hand-selected into a training set (n=18) and test set (n=4; B10, B10, T04, T18) [3]. The composition of the test set was changed, once the descriptors were selected for the final equations in order to get rid of any bias and compensate the small amount of compounds. The selection criterion for MLR equations was established with an R2 ≥ 0.4 [3].

Note, an additional fully automated 2D-QSAR approach (Ochem) is available online at https://ochem.eu/home/ show.do) [60]. In addition, the Supplementary Material presents another approach based on more physical grounds: the general regularization theory for ill-posed inverse problems was employed as an alternative to the MLR statistic technique in order to obtain a more general and robust linear QSAR equation. To this end, a smooth dependence from the activity data on the molecular descriptors was the only a priori assumption needed.

RESULTS

Synthesis

The organic synthesis resulted in yields in the range of milligrams and hereafter each compound could be crystallized and characterized. The melting points of all twenty synthesized compounds were each measured with a fully automated device known as EZ-Melt MPA120 and left uncorrected. Reaction monitoring was performed by thin layer chromatography (0. 2 mm pre-coated silica gel 60 F254 plates from E. Merck). 1H-NMR spectra were collected on a Varian Oxford (400 MHz) and 13C NMR (100 MHz) instruments. Chemical shifts were recorded in part per million (tetramethylsilane, δ =0) in deuterated dimethylsulfoxide; J values are represented in Hertz unities. Mass spectrometry analyses were carried out using a JEOL JMS-700 spectrometer under fast atom bombardment [FAB(+)]. Starting reagents were purchased from Sigma-Aldrich and the reagents did not undergo any purification before their use in the corresponding reaction procedures.

For more details concerning the spectrometric characterizations see the Supplementary Material.

Docked Ligand-PFOR Binding Analysis and 2D QSAR Modeling

Upon inspection of the docked ligand-enzyme poses, the interacting amino acids were identified (Table ) to answer questions about the possible activity enhancements related to chemical derivatization. This way we gained mechanistic insight and observed whether or not all ligands have one common binding mode. The proposed poses were compared with closely related crystallographic data to avoid a common pitfall, namely multiple binding modes, see Supplementary Material, Figs. , ) [3]. In theory, the NTZ derivatives could replace TPP in a competitive way. Actually, none showed higher binding energies than TPP itself. On the contrary, their affinities at best are found to be 1000-fold lower than TPP, so they cannot replace TPP as stronger binders to PFOR at all. Our computed finding is in excellent keeping with the postulated binding mechanism proposed by Hoffman et al. [32].

The 2D-QSAR model for biological response was build using stepwise forward multiple regression using all the calculated descriptors (as independent variables) and pIC50 values of biological response (as a dependent variable). Stepwise forward multiple regression adds descriptor one by one in the regression model until there is no significant improvement in training set R2 value. This analysis led to a multiple regression QSAR model with reasonable statistical parameters using one subset of four descriptors. Table reports the descriptors, regression coefficients and statistical parameters associated with the developed 2D-QSAR model (equation SA-1). The observed and predicted pIC50 values of training and test sets by using the 2D-QSAR model are shown in Fig. ().

The model revealed major role of BTZ_indicator descriptor i.e. 43.9% in determining Entamoeba histolytica activitiy, and it is inversely proportional to activity data. This suggests that compounds with thiazole scaffold are preferred over benzothiazole to attain higher potency against E. histolytica.

pIC50 = - 20.72 + 0.96 aasC_Cnt + 0.50 aOm_Cnt +

(equation SA-1)

3D-QSAR Equations for Activity Prediction

In principle, the QSAR equations predict new activities for hitherto unseen structures and therein support the decisions to take for a new cycle of synthesis, tests and improved design, although in this paper we focus more on the diagnostic strength of QSAR and possible shortcomings [3, 41, 42].

pIC50 = -3.85 + 2.17 pK (equation JL-1)

Analyzing the 3D-QSAR equations for the tilted compounds (equation JL-1, Fig. ) yields the following design message: the higher the descriptor values for a given compound – the more active it will become, since all three linear descriptors are positive with similar contribution weight. The first variable represents molecular acidity (pKa) which is formally located on the common exocyclic “>N-H” group and originated on the acidifying electron-withdrawing decoration of the mesomeric systems. The third one (lipole) is calculated in analogy to the dipole moment (DM) when molecular lipophilicity is broken down to all atoms of the molecule and assigned as fragmental atomic values to feed the dipole formula and thereby becoming what is called the lipole measure. The statistical analyses yield the following linear determination coefficients of R2Train=0.85, R2Test=0.96 with p<0.0001 (equation JL-1, Fig. ). In each case there is no risk of overfitting thanks to the 6:1 ratio between numbers of molecules versus descriptors [3]. The test set was hand-selected in the first place, but then during 5D-QSAR its composition was unattended, and the R2 did neither drop nor rise significantly (R2 ranged between 0.8 and 0.9, Fig. ).

In addition we present numerical results for the entire compound set, adding the four preselected compounds (B10, B30, T04, T18) of the test set. As can be noted for equation JL-1, the overall performance of our final QSAR model does by no means change, yet the values remain constant applying the LOO to the entire set of 22 compounds (see S).

Higher Dimensional QSAR

Concerning the higher dimensional QSAR (5D-QSAR), the model quality (3 equations JL-2, Fig. (), Bottom) is not significantly improved over the canonical QSAR linear equations (equation SA-1; R, n=18, R, n=4) or the 3D-QSAR equation (equation JL-1; R, n=18, R, n=4).

Inspecting the QSAR Models for Possible Pitfalls

We performed all QSAR studies on a training set (n=18) leaving four hand-selected molecules for the test set (B10, B30, T04, T18, Table ).

After inspection of two online resources (Marvin, SPARC [76, 77].) mayor inconsistencies were detected for the extended mesomeric systems concerning nitro groups [76] (pitfall: non constant “constants” [3]).

According to Table , seven out of twenty compounds possess higher activity than NTZ or TIZ in the in vitro susceptibility assays (T03, T04, T05, T06, T07, T08, and T09). The activity spreads over three orders of magnitude (pIC50 values) which complies with the rule of thumb concerning the data range pitfall [3]. However, the data set falls short of expectation concerning the covered chemical space, which is a serious setback when only a restricted variation of synthetic paths limit the functional enrichment of the produced compounds (cf. pitfall: data size and variety [3]). Despite obvious structural differences (cf. scaffold types of either thiazoles or benzothiazoles) we fused both clusters following the more pivotal constraint of dealing with a small data size and the limited chemical space therein (here: nTotal=22) [3].

In general, QSAR studies must be very carefully performed, i.e. under a balanced protocol (bias of training set, size and chemical diversity and activity range, representative internal and external test set, outlier handling) to ensure a meaningful interpretation of the chosen descriptors and require understanding of the mathematical operations that takes time and experience [78-84]. While practicing 2D- or 3D-QSAR modeling certain problems arose and were collected (Table ).

DISCUSSION

The application of computer aided drug design tools demonstrates the (more) rational approach to drug candidate development using medicinal chemistry concepts to support drug research. To guide new efforts in the near future into promising new chemical space, we predicted the theoretical activity of a new promising compound, see Supplementary Material (Fig. ).

Mechanistic insight comes from the docking studies indicating that NTZ and its derivatives do not replace cofactor TPP as a PFOR ligand at the binding site as proposed by Hoffman et al. [32].

As an asset, the present work sheds also light on the implications and limitations which the medicinal chemists and QSAR practitioners should be aware of as a sine qua non condition. If the physicochemical data space is huge, the linearity may be at stake (outliers, activity cliffs) whereas, if it becomes smaller and smaller or even ill-explored – like in the present case with only ether, nitro and halogen substituents – then QSAR prediction falls short of expectation without remedy. Apparently, very few data points are addressed and new structure predictions tend to become either intrapolations with no fresh ideas or extrapolations into uncovered chemical space.

Another challenge becomes the sheer number of descriptors available for QSAR model generation which makes it more probable that solutions are found “by chance”, which is obviously more the case in automated, unattended QSAR approaches with hundreds to thousands of descriptors [85]. In contrast, the meaningful descriptors can be hand-selected from the very beginning to build stepwise forward QSAR models, i.e. reiteratively increase the number of descriptors bit by bit (stepwise forward procedure). QSAR studies carried out without considering the active conformation compounds (tautomers, ionic forms, prodrugs) are omitting critical information reducing the quality of the input data needed for equations which in turn lose predictive power, ending up in an incorrect molecular design [3, 81, 82].

As an internal validation we used the final QSAR equations to predict the pIC50 values of the designed compound applying our linear equations, see Supplementary Material (Figs. ). Equation SA-1 yields 8.3 while equation JL-1 yields 8.4 which are practically the same.

CONCLUSION

The synthetic afford yielded seven new derivatives with a ten-fold or even stronger antiprotozoal activity than nitazoxanide. Concerning the relationships between structures and activities only mechanistically interpretable and fast to calculate descriptors were used in the final linear equations to predict even better candidates. Finally, we devised a hitherto unseen molecule with a 100-fold higher potency than NTZ, see Supplementary Material. The external validation was achieved upon manual and Random separation of molecules into a training and test set to guarantee representativeness despite the small input data size. This downside and other pitfalls were presented and solutions thereof proposed if possible. In addition, mechanistic insight came from the docking studies indicating that NTZ and its derivatives do not replace cofactor thiamine pyrophosphate as a pyruvate ferredoxin oxidoreductase binder at the active site as proposed in the scientific literature.

Table 1a

Listing of molecular structures of NTZ and its eleven thiazole derivatives (T series). The inhibition concentration pIC50 (against E. histolytica) is given and used as input data for QSAR. Note: tizoxanide (TIZ) is the hydrolysis product of NTZ: deacetyl-nitazoxanide. Note: recently, T17 without E. histolytica activity data was used in a published study (cf. Table 1 on page 1627 in [33]). And also recently T18 appeared without E. histolytica activity data (cf. Fig. 1 on page 3159 in another article by GNV [34]). Molecules of test set are denoted by asterisk (*).

Table 1b

Listing of molecular structures of NTZ and its nine benzothiazole derivatives (B series). The inhibition concentration pIC50 (against E. histolytica) are given and used as input data for QSAR. Note: B01 and B02 combined with E. histolytica activity data were already documented by the group leader GNV (cf. Table 1 on page 3169 in [27]. Molecules of test set are denoted by asterisk (*).

Table 2

Matrix of interactions representing Docking results (H-bonds, donors or acceptors). The gray boxes represent the interaction with the corresponding amino acid of PFOR model of E. histolytica. The experimental pIC50 and the computed ΔGbinding energies are included in the right side (T03 doesn’t show qualitative relationship between both values).

Table 3

Statistical parameters and descriptor definitions of the regression model (equation SA-1). Cross-validation standard error: cvSE; external-validation standard error: predSE; cross-validated Zscore: Zscore_cv; cross-validated alpha: alpha_cv.

Statistical Parameter	Value	Statistical Parameter	Value	Statistical Parameter	Value
Train/Test (n)	18/4	Q2F1	0.72	cvSE	0.97
Descriptors (k)	4	Q2F2	0.70	predSE	0.68
R2	0.80	Q2F3	0.69	Zscore	3.99
Q2 (LOO*)	0.55	r2m	0.39	Zscore_cv	2.42
pred R2	0.72			alpha	0.00
SEE	0.65	r2m (LOO*)	0.48	alpha_cv	0.01
Descriptor		Definition
aasC_Cnt		Count of atom-type E-State::C:-
aOm_Cnt		Count of atom-type E-State::O-
ssNH_Sum		Sum of atom-type E-State: -NH-
BTZ indicator		A binary variable either 0 or 1 to indicate the absence or presence of a benzothiazole scaffold in a given molecule

Note: (LOO*) is a particular parameter only applied to the initial training set (18 compounds).

Table 4

Listing of detected QSAR shortcomings and pitfalls which were reported in the literature [3].

Pitfalls	Comments
Small simple and limited chemical variability	It does not exist an ideal sample number size for QSAR, but is clear with larger sample size, the results become more representative. Here, the compound number is small (n=22), and the chemical variety is so limited. Basically it consist in ester (-O-CO-), ether (-O-), nitro (-NO2) and chloride, all of which occupy different positions of the phenyl ring of the molecules.
Composition of training and test sets	With a small sized series, the distribution of chemical items may not be representative in terms of the activity and chemical variability.
Meaningless descriptor selections	A common QSAR descriptor is pKa with high correlation to biological activity, which is not the case of dipole moment (DM, dipole). DM takes on different values with changing conformations. When DM is calculated for artificially held planar molecules, it is loaded onto the Z-axis only, while even small torsional changes are reflected by huge changes in DM values.
not constant coefficients and constants	The calculated pKa values are not equal to literature reports. It takes different values in different programs; moreover, each software consider different ionized forms, making difficult its selection/consideration for QSAR equations.Depending on the software, DM takes on different values due to normalization of input data.Certain descriptors like molar refractivity (MR) are very similar in magnitude in most programs. In contrast, polar surface area (PSA) should change according to the conformation of the molecule which implies that PSA is a 3D descriptors [44]. Other programs, however, calculate a “flat” PSA based on 2D data (atom connectivities and radii).
Starting geometries for 3D-QSAR	Albeit, the active conformation is not necessarily identical to the observed crystal structure, and since no NTZ-PFOR complex has been solved, two pieces of information were taken into account to assess the active conformation of the NTZ scaffold: (1) its crystallographic record deposited in CCDC [47] and (2) the final pose of NTZ docked into the ligand binding site of the cofactor TPP-PFOR complex. As a direct result, both geometries are practically the same, see Supplementary Material (Fig. SD2-C).
Errors of descriptor calculations (acidity, dissociation)	The experimental acidity value of NTZ is reported as pKa≈6 [32] for the conjugated acid / neutral thiazole system ([B-H+] / [B]) which corresponds to approx. 90% neutral species under physiological conditions. The calculated value, pKa≈8 [76], however, inverts the cationic/nonionic portions (10% neutral species). With no experimental value at hand, the (wrong) cationic forms would have been taken as input for the QSAR and docking studies.
Lipole-dipole collinearity	The algorithm of Lipole calculation is derived from the dipole moment equation (DM = q*r, q = atomic partial charge, r = VDW atomic radius) and atomic lipophilic values replace atomic partial charges. Despite different scale and units (charges and lipophilic fragments, same VDW radii), the equal calculation protocol generates collinearity.
Linearity hypothesis	The a priori assumption of linearity might be the main drawback in QSAR studies. Since data sampling is not complete, because no scientist would seek to explore the weaker, less active or more toxic data segments, it is often not clear if linearity is a first principle of nature or just appears due to insufficient data spread. Outliers and activity cliffs are first signs of nonlinear relationships between independent variables and response (biological activity, dependent variable).
Ligand based alignment (LBA)	The X-ray (crystal) conformation of NTZ may not constitute the biological active conformation. The hitherto unsolved structure of the NTZ-PFOR complex constitutes a disadvantage in case of higher dimensional QSAR where reliable conformational data is required. Results based on 2D descriptors (connectivities, drawings, SMILES, etc.) do not need special information while ligands can be superposed on their more rigid substructure or common scaffold (LBA).
Multiple solutions	We generated different equations based on different conformations and methodologies. It is not clear whether modeling based on NTZ X-ray conformation reflects realistic molecule geometries for binding site interaction, because the NTZ liganded binding site complex has not been elucidated. According to the pKa value of NTZ (here: 6.2, located on exothiazolic N amide), it can be inferred that all molecules treated here, present their activity at anionic form. Then, a new QSAR equation generation step based on descriptors calculated considering anionic compounds (without H at exothiazolic N amide, same training set), give us smaller R2 values that those obtained with X-ray data. The ideal case is considering the anionic form directly related with biological activity, because a small structural change can be reflected in huge descriptor magnitudes differences. This last QSAR equation generated with ionic compounds (pIC50= -2.36+2.28 pKa+0.17DM + 0.62 Lipole; R2Train=0.74, Q2=0.65, r2m=0.37, n=18; R2Test=0.75, Q2F1=0.75, Q2F2=0.75, Q2F3=0.75, r2m=0.68, n=4) could be seen as a poor predictive equation, but it becomes a better reflection of the biological behavior of our molecules.
Prodrugs and active metabolites	Some publications describe NTZ as a prodrug, albeit the biological activity of NTZ itself has been reported, too. Nevertheless, upon hydrolysis of the acetyl group, the metabolite TIZ shows comparable antiprotozoal potency.
Incompatible concepts and contradictions(chance correlation)	Sometimes, linear equations in 2D QSAR include conformation-dependent descriptors in a way where spatial information about structural requirements for the ligands and the binding site remains unknown. Hence, conformation-dependent descriptors contribute to establish the “rules” governing the relations between structures and activities, without any reason to be present in the equation except for chance correlations: “… because the relevant features only appear in molecules that also contain the wrong features” [3].