Structure-Based Design and Pharmacophore-Based Virtual Screening of Combinatorial Library of Triclosan Analogues Active against Enoyl-Acyl Carrier Protein Reductase of Plasmodium falciparum with Favourable ADME Profiles

Cost-effective therapy of neglected and tropical diseases such as malaria requires everlasting drug discovery efforts due to the rapidly emerging drug resistance of the plasmodium parasite. We have carried out computational design of new inhibitors of the enoyl-acyl carrier protein reductase (ENR) of Plasmodium falciparum (PfENR) using computer-aided combinatorial and pharmacophore-based molecular design. The Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) complexation QSAR model was developed for triclosan-based inhibitors (TCL) and a significant correlation was established between the calculated relative Gibbs free energies of complex formation (∆∆Gcom) between PfENR and TCL and the observed inhibitory potencies of the enzyme (IC50exp) for a training set of 20 known TCL analogues. Validation of the predictive power of the MM-PBSA QSAR model was carried out with the generation of 3D QSAR pharmacophore (PH4). We obtained a reasonable correlation between the relative Gibbs free energy of complex formation ∆∆Gcom and IC50exp values, which explained approximately 95% of the PfENR inhibition data: pIC50exp=−0.0544×∆∆Gcom+6.9336,R2=0.95. A similar agreement was established for the PH4 pharmacophore model of the PfENR inhibition (pIC50exp=0.9754×pIC50pre+0.1596, R2=0.98). Analysis of enzyme–inhibitor binding site interactions suggested suitable building blocks to be used in a virtual combinatorial library of 33,480 TCL analogues. Structural information derived from the complexation model and the PH4 pharmacophore guided us through in silico screening of the virtual combinatorial library of TCL analogues to finally identify potential new TCL inhibitors effective at low nanomolar concentrations. Virtual screening of the library by PfENR-PH4 led to a predicted IC50pre value for the best inhibitor candidate as low as 1.9 nM. Finally, the stability of PfENR-TCLx complexes and the flexibility of the active conformation of the inhibitor for selected top-ranking TCL analogues were checked with the help of molecular dynamics. This computational study resulted in a set of proposed new potent inhibitors with predicted antimalarial effects and favourable pharmacokinetic profiles that act on a novel pharmacological target, PfENR.


Introduction
Malaria is a major disease in tropical areas, where about 40% of the world's population lives. This disease is responsible for about 0.6 million deaths per year, 98% of which are in Africa. According to the WHO, children under 5 years of age account for 77% of all malaria deaths [1]. Plasmodium falciparum (Pf ) infection is the most virulent form of malaria, accounting for most malaria cases worldwide [1]. The significant spread of resistance to traditional antimalarial drugs underscores the need to develop new potent and orally bioavailable antimalarial agents, which are also active against existing drug-resistant P. falciparum strains [2].
Recently reported cases of resistance of Pf to artemisinin and combined therapies (ACT) in Rwanda, Uganda, and Eritrea [3,4] as well as in Burkina Faso and Mali [5] worsen the picture and make the discovery of new antimalarials for Africa more urgent. Finally, almost 400 years after the discovery of quinine (Cinchona bark) originating in South America and 40 years after the introduction of artemisinin (Artemisia annua) from Asia, eyes are fixed on Africa and finding a novel molecule from the rich pharmacopoeia of this continent, which will help overcome the problem of the resistance of Pf to artemisinin and its derivatives. Currently, numerous validated pharmacological targets for Pf inhibition have been identified at different stages of the complex parasite life cycle, many of them not yet involved in the parasite's resistance [6,7]. In fact, the treatment of parasite infection in the erythrocytic stage is preferred for most approved antimalarials, in contrast to the liver stage, where administered medical preparations serve primarily for prophylactic purposes [8,9].
The synthesis of fatty acids is an essential process for all living organisms. Fatty acids play a critical role in the supply of metabolic precursors of biological membranes and represent an important form of metabolic energy production. Attention is being given to inhibition of the type II fatty acid biosynthesis pathway (FAS-II), which is different from the related FAS-I pathway present in mammalian cells. FAS-II of Pf appears to hold significant promise in the discovery of novel drugs, which are expected to be less toxic to humans due to the difference between these two distinct pathways [10][11][12][13]. Enoyl-acyl carrier protein reductase (ENR) is a crucial enzyme in the type II fatty acid synthesis pathway of many pathogens, including Plasmodium falciparum, the etiologic agent of malaria. The final reaction in the FAS-II pathway of Pf is catalysed by enoyl-acyl carrier protein reductase (Pf ENR), which mediates the NADH-dependent reduction of the trans-2-enoyl-acyl carrier protein [14,15].
Triclosan (TCL1, Figure 1), a common antibacterial agent incorporated into a variety of consumer products, is an uncompetitive inhibitor of purified Pf ENR, which has demonstrated inhibitory potency against P. falciparum parasites cultured in vitro IC exp 50 = 73nM [16][17][18]. The TCL1 derivatives ( Figure 1) are non-toxic to mammalian cells [17,19] which do not contain ENR homologues. The human exposure routes, metabolism, and toxicity of triclosan have recently been reviewed with a conclusion that adverse biological effects of triclosan and its degradation products need to be investigated [20]. Freundlich previously reported that TCL1 analogues inhibit both Pf ENR enzymes and Mycobacterium tuberculosis ENR enzymes [21][22][23][24][25]. Pf ENR is a well-known validated pharmacological target of antimalarials [26]. The existence of X-ray crystal structures of TCL1 and its analogues plus the NAD + cofactor bound to Pf ENR facilitate the possibility of computer-assisted drug design of antimalarial agents, which target lipid synthesis in the plasmodium parasite [18,21]. In the bound conformation, the aromatic rings of TCL1 are orthogonal to each other, the phenol ring mimicking the enolate intermediate of the reaction catalysed by the ENR. Simultaneously the TCL1 is face to face π-stacked with NAD+, T-shaped π-stacked with the catalytic Tyr 277, and forms hydrogen bonds to the cofactor and the tyrosine, respectively [21]. Replacements of chlorine atoms at positions 5, 2 and 4 on the TCL1 diarylether scaffold by larger groups resulted in energetically favourable interactions with both the enzyme and the co-factor [22,23,27,28]. Synthesis of 4 -derivatives of TCL1 caused minor gains in the inhibitory potencies towards Pf ENR [27], while the analogues substituted in the 2 position reached only the micromolar range of the IC 50 s [23]. Sullivan et al. [28] described the synthesis and inhibitory activity tests of 5-alkyl-substituted TCL1 derivatives toward purified Pf ENR. Chibber et al. [22] studied TCL analogues with mainly hydrophilic 5-substituents in carbon 5, which were less potent than TCL1. Freundlich et al. [21] reported on the synthesis and activity tests of TCL1 analogues containing hydrophobic substituents at position 5 intended to better 'root' the inhibitor at the active site of Pf ENR. Their results demonstrated the possibility of improving the in vitro antimalarial potency of TCL1 by replacing the suboptimal 5-chloro group with larger hydrophobic moieties. The TCL1 scaffold has inspired numerous initiatives to design potent Pf ENR inhibitors in a search for new antimalarials [29][30][31]. In our previous work, a QSAR model of TCL binding to PfENR was built using the LUDI scoring function to allow the in silico screening of a virtual combinatorial library of 120 TCL1 analogues with favourable predicted ADME profiles [32]. The best designed TCLs reached predicted IC pre 50 slightly over 10 nM (Figure 1). So far, few 3D-QSAR pharmacophore (PH4) studies have been reported for the inhibition of PfENR by TCL. The PH4 pharmacophore was used for TCL1 in complex with Pf ENR and based on receptor coordinates recovered from the X-ray structure 3F4B (Protein Data Bank) [33]. Unfortunately, this PH4 model is spatially limited, since 5 -Cl is not large enough to fill the surrounding hydrophobic pocket of the enzyme active site. The search for a reliable PH4 model taking into account the reported SAR of Pf ENR inhibition for the known TCL inhibitors [9,[14][15][16]18] seems tricky, as illustrated in the following comparative analysis of the role of 5 -substituents: TCL1 (R 1 , R 2 , R 3 ≡ Cl, Cl, Cl; 73 nM), TCL10 (R 1 , R 2 , R 3 ≡ CH 2 -Ph, Cl, Cl; 71 nM), TCL11 (R 1 , R 2 , R 3 ≡ (CH 2 ) 2 -Ph, Cl, Cl; 76 nM), and TCL15 (R 1 , R 2 , R 3 ≡pF-Ph, Cl, Cl; 38 nM). Extension of the -(CH 2 ) n -linker to fill the hydrophobic pocket next to C-5 is not advantageous, while the p-fluorophenyl of TCL15 increases the potency twice, in contrast to linkers. Therefore, the design of new more potent TCLs against Pf ENR requires careful application of the 3D-QSAR PH4 model, attention to the active conformation, and also exploration of the pockets of the active site for the positions C2 and C4 of the TCL scaffold.  Table 1) and 2D interaction diagram with residues from the active site of PfENR (PDB entry 2OOS) [21]. Atom colouring: C-grey, O-red, Cl-green. (C) triclosan analogues that were previously predicted to be potent inhibitors of PfENR [6,32]. Table 1. Training and validation sets of TCL inhibitors used for the preparation of QSAR models of PfENR inhibition. The experimental activities were taken from [21,23,25,27].
In this work, we further develop predicted inhibitory potencies of TCLs against Pf ENR by extended combinatorial design and a virtual screening approach: -We built and validated a Hansch-type QSAR model to correlate the relative Gibbs free energy (rGFE) of TCL binding to Pf ENR with the observed inhibitory potencies IC exp 50 (identification of the bound conformation of TCL inhibitors). - The robustness of the built model is confirmed by a 3D-QSAR pharmacophore (PH4) model based on the bound conformations of the training set of inhibitors (generation of the PH4 model).
-In addition, we built a virtual library of 33,480 analogues of TCLs and screened them with the help of PH4 to identify hits for inhibition of Pf ENR (virtual screening). -Finally, the ADME profile of the best designed analogues was predicted and compared with those of current antimalarial drugs and compounds undergoing clinical trials (postprocessing step 1). -Each of the hits was cross-checked for predicted inhibitory potency by computed rGFE of the formation of the Pf ENR-TCLx complex, which led to the identification of a handful of prospective novel TCLs that exceeded the potency of known TCLs against Pf ENR (postprocessing step 2). - The top 5 TCL hits and TCL11 underwent molecular dynamics simulations to explore the stability of the Pf ENR-TCLx complexes (postprocessing step 3).  Table 1. Training and validation sets of TCL inhibitors used for the preparation of QSAR models of PfENR inhibition. The experimental activities were taken from [21,23,25,27].

Training and Validation Sets of TCL Inhibitors
The training set (20 molecules) and validation set (4 molecules) of triclosan analogues that inhibit Pf ENR of the malaria parasite (Table 1) used in this study were selected from the literature [21,23,25,27]. The inhibitory potencies of these derivatives cover a sufficiently broad range of activities to allow a reliable QSAR model to be built 38 nM ≤ IC exp 50 ≤ 10 5 nM .  (Table 1), an enzyme-inhibitor complex Pf ENR-TCLx was constructed by in situ modifications of the refined template crystal structure of the complex Pf ENR-TCL11 (PDB entry code 2OOS [21]) as described in Section 4. Furthermore, the relative Gibbs free energy (rGFE) of the complex formation (∆∆G com ) and its components was calculated for each of the 24 geometry-optimised Pf ENR-TCLx complexes for the triclosan derivatives ( Table 2; see footnote of Table 2 for definitions of the quantities). A QSAR model that correlates the computed rGFE with the experimental inhibitory potencies of TCLs (pIC exp 50 = −log 10 IC exp 50 [17]) explained 95% of the variation in the observed activities. The resulting linear regression equation is shown in Table 3, Equation (2). Relatively high values of the regression coefficient R 2 , leave-one-out cross-validated regression coefficient R 2 xv and Fischer F-test of the correlation suggest a strong relationship between the 3D model of inhibitor binding and the observed inhibitory potencies of TCL [21]. Therefore, structural information derived from the 3D models of Pf ENR-TCLx complexes can be expected to lead to a reliable prediction of the inhibitory potencies of Pf ENR for new TCL analogues based on QSAR Equation (2) ( Table 3). To gain a better understanding of the binding affinity of TCLs to Pf ENR, we have analysed the enthalpy of complexation in gas phase ∆∆H MM by correlating it with the pIC exp 50 . The validity of the resulting linear correlation (for statistical data from the regression, see Table 3, Equation (1)) allowed the assessment of the contribution of inhibitor-enzyme interactions (∆∆H MM ) when the solvent effect and loss of inhibitor entropy upon binding to the enzyme were neglected. This correlation explained about 94% of the pIC exp 50 data variation and underlined the role of interatomic interactions in determining the binding affinity of the ligands. These statistically significant correlations allowed us to reveal the bound conformations of the TCLs at the Pf ENR binding site and allowed the definition of the PH4 pharmacophore.
The statistical data confirmed the validity of the correlation Equations (1) and (2) shown in Figure 2. The ratio pIC  Table 3). Thus, the regression Equation (2) and computed ∆∆G com (Table 2) can be used to predict inhibitory potencies against Pf ENR for new TCL analogues, provided that they share the same binding mode as the training set of triclosan derivatives TCL1-20.  [22,24,26,28]. g The ratio of predicted and experimental half-maximal inhibitory concentrations pIC pre 50 /pIC exp 50 . pIC pre 50 = −log 10 IC pre 50 was predicted from the computed ∆∆G com using the regression Equation (2) for Pf ENR shown in Table 3.    Table 3) of the training and validation sets [21]. The data points from the validation set are shown in red.

Binding Mode of TCLs to PfENR
The binding mode of TCLs derived from the crystal structure [21] and the QSAR model are illustrated in Figures 3 and 4. The interactions of the TCL inhibitor at the active site of PfENR present in the reference X-ray structure of the TCL11 complex (Figure 3, [21]), namely the -stacking between the phenol moiety and the nicotinamide ring of ∆∆ ∆∆   Table 3) of the training and validation sets [21]. The data points from the validation set are shown in red.

Binding Mode of TCLs to Pf ENR
The binding mode of TCLs derived from the crystal structure [21] and the QSAR model are illustrated in Figures 3 and 4. The interactions of the TCL inhibitor at the active site of Pf ENR present in the reference X-ray structure of the TCL11 complex ( Figure 3, [21]), namely the π-π stacking between the phenol moiety and the nicotinamide ring of NAD + , are conserved. Furthermore, hydrogen bonds (HB) involving the hydroxyl group of TCL and Tyr277, as well as the ribose ring of NAD + (Figures 3 and 4), are also preserved. One of the main contributions of TCL derivatives to the enzyme-inhibitor interaction in PfENR involves the hydrophobic pocket formed by the residues Tyr267, Val274, Tyr277, Tyr375, Phe368, Ala269, Ala320, Ala372, and Ile369 [21,32]. As we can see in Figure 4, for inhibitor TCL15 the hydrophobic R substituent Ph pF reaches the hydrophobic pocket and contributes significantly to the stabilising interactions at the active site of PfENR. The R group of TCL analogues is surrounded by active site residues Asn218, Val222, Met281, and substrate binding loop residues Ala319 and Ile323.   [21]). The carbon atoms are coloured yellow; the residue side chains and carbon atoms of the NAD cofactor are coloured cyan and purple, respectively. Interaction colour code: hydrogen bonds (green), hydrophobic contacts of alkyl groups (pink), π-π interactions (dark pink). (Right) Molecular surface of the active site of PfENR. Surface colouring legend: red, hydrophobic; blue, hydrophilic; white, intermediate.

Further enrichment of enzyme-inhibitor interactions by new substitutions in TCL1
analogues with the same mode of PfENR binding will lead to improved potencies of the new structurally similar compounds, provided that the predictive power of correlation Equation (2) also extends beyond the range of activities of the training set. The predictive ability of our QSAR model relying on the computed ∆∆G was also confirmed by the PH4 pharmacophore model of the PfENR inhibitory activity.

3D-QSAR Pharmacophore Model of PfENR Inhibition
The virtual screening process is operated with a 'chemical space navigator' (Pharmacophore, PH4) built from the bound conformations of the most active of 20 TCLs recov-  [21]). The carbon atoms are coloured yellow; the residue side chains and carbon atoms of the NAD cofactor are coloured cyan and purple, respectively. Interaction colour code: hydrogen bonds (green), hydrophobic contacts of alkyl groups (pink), π -π interactions (dark pink). (Right) Molecular surface of the active site of Pf ENR. Surface colouring legend: red, hydrophobic; blue, hydrophilic; white, intermediate.
One of the main contributions of TCL derivatives to the enzyme-inhibitor interaction in Pf ENR involves the hydrophobic pocket formed by the residues Tyr267, Val274, Tyr277, Tyr375, Phe368, Ala269, Ala320, Ala372, and Ile369 [21,32]. As we can see in Figure 4, for inhibitor TCL15 the hydrophobic R 1 substituent Ph −pF reaches the hydrophobic pocket and contributes significantly to the stabilising interactions at the active site of Pf ENR. The R 2 group of TCL analogues is surrounded by active site residues Asn218, Val222, Met281, and substrate binding loop residues Ala319 and Ile323.
Further enrichment of enzyme-inhibitor interactions by new substitutions in TCL1 analogues with the same mode of Pf ENR binding will lead to improved potencies of the new structurally similar compounds, provided that the predictive power of correlation Equation (2) also extends beyond the range of activities of the training set. The predictive ability of our QSAR model relying on the computed ∆∆G com was also confirmed by the PH4 pharmacophore model of the Pf ENR inhibitory activity.

3D-QSAR Pharmacophore Model of PfENR Inhibition
The virtual screening process is operated with a 'chemical space navigator' (Pharmacophore, PH4) built from the bound conformations of the most active of 20 TCLs recovered from the Pf ENR-TCL1-20 complexes that served to develop the QSAR model. Unlike direct pharmacophore generation from small free (not bound) molecules, the HypoGen algorithm in the 3D-QSAR pharmacophore protocol of Discovery Studio (DS) [34] generates a PH4 in order to estimate the biological activity (IC 50 ) of the 20 training set compounds. Their IC exp 50 values span 3.42 orders of magnitude (from 38 to 10 5 nM), which is very close to the recommended 3.5 for reliable PH4 generation. To generate PH4, the constructive step of the algorithm selects compounds with IC exp 50 ≤ 1.25 × 38 nM in order to build the starting PH4 features; only TCL15 (IC exp 50 = 38 nM) is retained for this purpose. The 3.5 orders of magnitude criterion also served to identify the TCL compounds to be considered inactive in the subtractive step of the PH4 generation (IC exp 50 > 38 × 10 3.5 nM = 120167 nM). Therefore, none of the other compounds from TCL1-20 were subtracted, which means that all the PH4 features generated from TCL15 were conserved for the processing of the hypotheses. Finally, during the last PH4 generation step, namely the optimisation phase, the score of the PH4 hypotheses is improved. According to a detailed description in the literature [35] as well as in Section 4, four main features were identified: hydrophobic (HYD) complementary to enzyme lipid pockets, hydrogen bond donor (HBD), and acceptor (HBA) binding to active site residues by acceptor and donor hydrogen bonds, respectively. Virtually screened molecules were mapped to these four pharmacophore features. The output of the HypoGen process is the top 10 PH4 hypotheses ranked according to their statistical data as presented in Table 4. Among them, the hypothesis cost is one of the most pertinent: the metrics defines the hypothesis total cost as a linear combination of three components, namely the weight, the error, and the configuration cost (see Section 4). It spans a range from 409.3 (Hypo1) to 1231.8 (Hypo10), all in a monotonic increasing order. Since the total cost is 803.5 (almost twice that of Hypo1), Hypo1 is the closest to the fixed cost and is very far from the other nine models. For the entire set of 10 hypotheses, two extreme additional costs are computed: the fixed cost (lower limit) of 15.2 and the null cost (upper bound) of 5423.2; the difference ∆ = 5408.1 between them is the first criterion of the analysis. The larger the value of ∆ (∆ > 70), the more meaningful are the top hypothesis models with a statistical significance greater than 90%. The second indicator is the difference between the null cost and the total cost of each hypothesis which is higher than 4190, indicating that all the 10 hypothesis models are very far from the null cost. The third indicator is the configuration cost of 9.55, far below 17, which is the value it should not exceed [36]. The next descriptor is the root-mean-square deviation (RMSD) of each hypothesis. The RMSD value of Hypo1 is equal to 6.6, the lowest. Hypo2 reached a RMSD of 9.5. Finally, the correlation coefficient (R 2 ) of the training set reached a value of 0.96 for Hypo1 and 0.93 for the following Hypo2. All generated hypotheses display HYD, HBA, and HBD features. Finally, the first hypothesis (Hypo1) was retained for the in silico screening of the virtual library of TCL analogues. The additional statistical parameters of the regression for pIC exp 50 vs. pIC  Figure 5 and is very close to a slope (1.02) of one and an intercept of 0.13 near the origin. The parameters reflect the guidance of the Organisation of Economic Cooperation and Development (OECD) about QSAR [37]. An additional verification is the mapping of the most potent inhibitor TCL15 to the Hypo1 features as shown in Figure 5. In the same figure, the detailed geometry and position of the Hypo1 PH4 features of Pf ENR inhibition are presented. PH4 was externally evaluated by a validation set of 4 derivatives, TCL21-24, that belong to the same interval of inhibitory potencies as the compounds of the training set.    For each of the 10 hypotheses, a confidence level of 98% is retained. The number of scrambled runs required is obtained using the following equation

Virtual Library Generation and In Silico Screening of PfENR Inhibitors
where the number of hypotheses (X) the total cost of which is below that of Hypothesis 1 (Hypo1), is equal to 0; the total number of HypoGen runs (initial + random runs) should be equal to 1 + 49.

Replacement of these values in the equation results in
The validation of the statistical significance of each hypothesis model is carried out using Fisher's randomisation method. The process consists of the randomisation of the observed biological IC exp 50 for the 20 training set compounds. It checks whether the PH4 model is not the result of an incidental correlation between IC exp 50 and the conformations of TCL compounds (CatScramble programme). As indicated above, for the confidence level of S = 98%, 49 sets are generated for 10 pharmacophore hypotheses from spreadsheets built by random relationships between the IC exp 50 and 3D structures of the TCL1-20. From the lowest cost among the 49 hypotheses for each of the 10 PH4 (Hypo1-Hypo10) listed in Table 4 (closest random), none shows a cost lower than the 409.3 of Hypo1. Randomisation confirms that Hypo1 is not produced by chance and is a true correlation supported by a statistically robust pharmacophore model for inhibition of Pf ENR by TCL with a confidence level of 98%. It shows predictive power as the QSAR complexation model on the ligand affinity descriptor ∆∆G com .
Next, we have carried out the computational design and selection of new TCL analogues with predicted elevated inhibition potencies against the malaria parasite Pf ENR, based mainly on the presence of the hydrophobic feature (HYD, Figure 5) in the best PH4 model at the substitution position of the R 1 group (see Table 1) and two hydrophobic aliphatic features located at the R 2 and R 3 positions. (HYD_AL-2 and HYD_AL-1, Figures 1 and 5).

Virtual Library Generation and In Silico Screening of PfENR Inhibitors
In silico screening of a virtual combinatorial library can lead to hit identification, as was shown in our previous work on inhibitor design [32].

Virtual Library
In a previous study on the antimalarial effects of TCLs [32], we have proposed a highly focused virtual combinatorial subset. Therefore, 40 smaller polar and larger aliphatic and aromatic hydrophobic R-groups (fragments, building blocks, reagents) were selected as suitable diverse building blocks for the position R 1 . Out of them, 23 and 11 smaller aliphatic reagents were chosen as building blocks for the positions R 2 and R 3 , respectively, which are known to be able to accommodate smaller substituents.
In this work, the virtual combinatorial subsets were extended to 60, 31 and 18 R-groups. The R-groups listed in Table 5 were attached to the positions R 1 , R 2 and R 3 of the TCL scaffold to form a combinatorial library of the size: R 1 × R 2 × R 3 = 60 × 31 × 18 = 33, 480 TCL analogues. This initial diversity library was generated from building blocks (chemicals) listed in the database of available chemicals [38].

In Silico Screening
The generated virtual diversity library of 33,480 TCL analogues was further screened for molecular structures that matched the 3D-QSAR pharmacophore models of Pf ENR inhibition to identify compounds that are likely to bind to this target. During virtual screening, 1000 conformers were generated for each of 33,480 TCL analogues. After screening, 219 TCLs mapped to at least 2 pharmacophoric features of the 3D-QSAR pharmacophore model Hypo1 ( Figure 5A,B) of inhibition of Pf ENR. These fitting analogues (PH4 hits) were then further evaluated by the complexation model. The computed relative GFE of the formation of the Pf ENR-TCLx complex, its components, and the predicted half-maximal inhibitory concentrations IC pre 50 calculated from the correlation Equation (2) ( Table 3) are listed in Table 6 (for a complete table, see Table S1, Supplementary Material).

Novel TCLs against PfENR
The results of the in silico screening of the virtual library of novel TCL analogues were also analysed in terms of the frequency with which each individual fragment appears in the R groups of the 219 virtual hits that match the PH4 model of inhibition of Pf ENR. The R 1 moiety of TCL that fills the hydrophobic pocket corresponding to the HYD pharmacophoric feature is a key determinant of improving the potency of TCLs against the Pf ENR. Five  1(18) contain chiefly small primary, secondary, and tertiary amines ( Figure 6). Table 6. Relative GFE calculated and its components for 20 best analogues of TCLs identified by in silico screening of the diversity combinatorial library. The numbering of the analogues concatenates the index of each substituent R 1 − R 2 − R 3 with the substituent numbers taken from  Table 2); b ∆∆G sol is the relative solvation contribution to ∆∆G com ; c ∆∆TS vib is the relative entropic (vibrational) contribution to ∆∆G com ; d ∆∆G com is the relative GFE change of the Pf ENR-TCLx complex formation ∆∆G com = ∆∆H MM + ∆∆G sol − ∆∆TS vib ; e pIC pre 50 is the logarithm of predicted inhibition potency towards Pf ENR calculated from ∆∆G com using correlation Equation (2), Table 3 (Tables 5 and 6). They include the following fragments at R 1 -position: 33: ···CH 2 -cyclohexyl; 58: ···n-hexyl; 59: ···n-heptyl; 60: ···n-octyl. At R 2 -position: 1: ···H; 8: ···NH 3 + ; 16: ···trimethylamine; 23: ···Ph; and at R 3 -position 1: ···H; 5: ···OH; 8: ···NH 3 + ; 9: ···C 2 H 5 . Due to the amino acid composition of the larger hydrophobic pocket, the R 1 -groups display preferences for bulkier nonpolar building blocks of quite variable topology. The R 2 -groups do not show a common character either, while the R 3 -groups show preference mainly for small polar or cationic substituents. In the specific case of Pf ENR, the hydrophobic contacts are expected to root aryl-substituted inhibitors at the active site. The presence of small polar and cationic substituents in the R 3 -position resulted in strong electrostatic attraction with the proximate phosphate groups of the cofactor nicotinamide adenine dinucleotide. The cationic substituents in the partially solvent-exposed R 3 -position caused stabilization of the bound TCL analogues by the effect of hydration. These substitutions led to an overall predicted increase in the binding affinity of the new TCL analogues to Pf ENR. The best TCL analogues TCL-60-08-08 and TCL-33-08-05 (Figure 7) show predicted half-maximum inhibitory concentrations of (IC    The role of hydrophobic residues of the active site pocket surrounding the R1-groups of TCLs is a key determinant of the predicted increase in the affinity of virtual hits. To gain a deeper understanding of enzyme-inhibitor interactions, the calculated interaction energy (E ) of the most potent training set residues with the PfENR was analysed in The role of hydrophobic residues of the active site pocket surrounding the R 1 -groups of TCLs is a key determinant of the predicted increase in the affinity of virtual hits. To gain a deeper understanding of enzyme-inhibitor interactions, the calculated interaction energy (E int ) of the most potent training set residues with the Pf ENR was analysed in terms of contributions from individual residues of the active site (Figure 8). To compare E int contributions of important active site residues to the total enzyme-inhibitor interaction energy for TCL1 (IC

Pharmacokinetic Profile of Novel TCL Analogues
Properties related to ADME (absorption, distribution, metabolism, and excretion) were estimated for virtual hits, as well as for some selected reference antimalarial drugs. The best designed TCL derivatives with low oral bioavailability due to poor water solubility and rapid phase II metabolism must be disregarded for possible use as antimalarials. All of the ADME-related properties shown in Table 7, such as the octanol-water partition coefficient; aqueous solubility; blood-brain partition coefficient; Caco-2 cell permeability; serum protein binding; number of likely metabolic reactions; and another 18 descriptors of the new analogues, were computed by the QikProp program [39], based on the method of Jorgensen [40,41]. Experimental data from more than 710 compounds including about 500 drugs and related heterocycles were used to produce regression equations that correlated experimental and computed descriptors, resulting in a fast and accurate prediction of the pharmacokinetic properties of drug-like molecules. Consistent with the unfavourable pharmacokinetic characteristics observed for the best active triclosan derivatives, the predicted oral bioavailability of the new TCL analogues ranges from 86% to 100%. Since a value greater than 80% is considered good, six analogues among the best 10 hits show a human oral absorption in the gastrointestinal tract (HOA) of 100%. Drug-likeness (#stars)-the number of property descriptors that fall outside the optimal value range determined for 95% of known drugs out of 24 selected descriptors computed by QikPropwas used as an additional general criterion for selecting compounds with good ADMErelated properties. The values of the best active designed TCLs against PfENR are compared with those computed for drugs used for the treatment of malaria alone or, in the case of artemisinin, drugs used in a combined therapy (ACT), or drugs currently undergoing clinical trials (Table 7). This comparison revealed that triclosan analogues are endowed with ADME-related properties which are comparable to most of the considered antimalarial drugs.

Pharmacokinetic Profile of Novel TCL Analogues
Properties related to ADME (absorption, distribution, metabolism, and excretion) were estimated for virtual hits, as well as for some selected reference antimalarial drugs. The best designed TCL derivatives with low oral bioavailability due to poor water solubility and rapid phase II metabolism must be disregarded for possible use as antimalarials. All of the ADME-related properties shown in Table 7, such as the octanol-water partition coefficient; aqueous solubility; blood-brain partition coefficient; Caco-2 cell permeability; serum protein binding; number of likely metabolic reactions; and another 18 descriptors of the new analogues, were computed by the QikProp program [39], based on the method of Jorgensen [40,41]. Experimental data from more than 710 compounds including about 500 drugs and related heterocycles were used to produce regression equations that correlated experimental and computed descriptors, resulting in a fast and accurate prediction of the pharmacokinetic properties of drug-like molecules. Consistent with the unfavourable pharmacokinetic characteristics observed for the best active triclosan derivatives, the predicted oral bioavailability of the new TCL analogues ranges from 86% to 100%. Since a value greater than 80% is considered good, six analogues among the best 10 hits show a human oral absorption in the gastrointestinal tract (HOA) of 100%. Drug-likeness (#stars)-the number of property descriptors that fall outside the optimal value range determined for 95% of known drugs out of 24 selected descriptors computed by QikProp-was used as an additional general criterion for selecting compounds with good ADME-related properties. The values of the best active designed TCLs against Pf ENR are compared with those computed for drugs used for the treatment of malaria alone or, in the case of artemisinin, drugs used in a combined therapy (ACT), or drugs currently undergoing clinical trials (Table 7). This comparison revealed that triclosan analogues are endowed with ADME-related properties which are comparable to most of the considered antimalarial drugs.   ); e hydrophobic portion of the solvent-accessible molecular surface, in (Å 2 ) (range for 95% of drugs: 0-750 Å 2 ); f total volume of the molecule enclosed by solvent-accessible molecular surface, in (Å 3 ) (range for 95% of drugs: 500-2000 Å 3 ); g number of rotatable non-trivial (not CX3), non-hindered (not alkene, amide, small ring) bonds (range for 95% of drugs: 0-15); h estimated number of hydrogen bonds that would be donated by the solute to water molecules in an aqueous solution. Values are averages taken over a number of configurations, so they can assume non-integer values (range for 95% of drugs: 0.0-6.0); i estimated number of hydrogen bonds that would be accepted by the solute from water molecules in an aqueous solution (range for 95% of drugs: 2.0-20.0); j logarithm of the partition coefficient between the n-octanol and water phases (range for 95% of drugs: −2 to 6.5); k logarithm of predicted aqueous solubility, logS. S in (mol·dm -3 ) is the concentration of the solute in a saturated solution that is in equilibrium with the crystalline solid (range for 95% of drugs: −6.0 to 0.5); l logarithm of the predicted binding constant to human serum albumin (range for 95% of drugs: −1.5 to 1.5); m logarithm of the predicted brain/blood partition coefficient (range for 95% of drugs: −3.0 to 1.2); n predicted apparent Caco-2 cell membrane permeability on the Boehringer-Ingelheim scale in [nm·s −1 ] (range for 95% of drugs: <25 poor, >500 nm·s −1 great); o number of likely metabolic reactions (range for 95% of drugs: 1-8); p predicted inhibition constants IC pre 50 of designed TCLs vs. Pf ENR.; q human oral absorption (1 = low, 2 = medium, 3 = high); r percentage of human oral absorption in the gastrointestinal tract (<25% = poor, >80% = high); * star in any column indicates that the property descriptor value of the compound falls outside the range of values for 95% of known drugs.

Molecular Dynamics Simulations
We also performed molecular dynamics simulations to check the stability of Pf ENRinhibitor complexes and the flexibility of the active conformations of the native ligand TCL11 and five of the best-designed TCL analogues at the active site of Pf ENR. Thus, the complexes obtained from the modification in situ of the reference inhibitor TCL11 and subsequent refinement through molecular mechanics were used as starting geometries for the molecular dynamics calculations using the Desmond software [42] (see Section 4). Table 8 presents the averages of total energy E tot and potential energy E pot of the systems during the 200 ns simulation. Figure 9 illustrates the simulated system and Figure 10 shows the time evolution of bound inhibitor properties, such as the root mean square deviation (RMSD) from the initial conformation, the radius of gyration (rGyr), number of intramolecular hydrogen bonds (intraHB) molecular surface area (molSA), solvent-accessible surface area (SASA) and the polar surface area (PSA). Protein-ligand interactions were explored throughout the simulation trajectory to identify specific interactions maintained during computation ( Figure 11). Finally, the interactions that occur more than 20.0% of the simulation time in the trajectory (0-200 ns) are shown on a detailed 2D schematic representation ( Figure 12). Moreover, we have superimposed the conformations of the ligands obtained following the minimisation of the complexes from molecular dynamics on those obtained by in situ modification of the TCL11 and MM refinement. Figure 13 illustrates these superimpositions and the respective RMSDs. The RMSD, rGyr, SASA, and the remaining parameters document that the pose of the bound new TCLs in the complexes with Pf ENR remains stable within the MD simulation time (Figure 10). The higher contribution of hydrogen bonding interactions of the new TCL analogues with the active site residues as compared to the predominantly hydrophobic interaction of the TCL11 reference compound suggests an elevated specificity of these virtual hits for the Pf ENR target (Figure 11).

Molecular Dynamics Simulations
We also performed molecular dynamics simulations to check the stability of inhibitor complexes and the flexibility of the active conformations of the native TCL11 and five of the best-designed TCL analogues at the active site of PfENR. Th complexes obtained from the modification in situ of the reference inhibitor TCL subsequent refinement through molecular mechanics were used as starting geomet the molecular dynamics calculations using the Desmond software [42] (see Section ble 8 presents the averages of total energy 〈E 〉 and potential energy 〈E 〉 of t tems during the 200 ns simulation. Figure 9 illustrates the simulated system and Fi shows the time evolution of bound inhibitor properties, such as the root mean squ viation (RMSD) from the initial conformation, the radius of gyration (rGyr), num intramolecular hydrogen bonds (intraHB) molecular surface area (molSA), solven sible surface area (SASA) and the polar surface area (PSA). Protein-ligand inter were explored throughout the simulation trajectory to identify specific interaction tained during computation ( Figure 11). Finally, the interactions that occur mo 20.0% of the simulation time in the trajectory (0-200 ns) are shown on a detailed 2 matic representation ( Figure 12). Moreover, we have superimposed the conforma the ligands obtained following the minimisation of the complexes from molecular d ics on those obtained by in situ modification of the TCL11 and MM refinement. Fi illustrates these superimpositions and the respective RMSDs. The RMSD, rGyr, and the remaining parameters document that the pose of the bound new TCLs in th plexes with PfENR remains stable within the MD simulation time ( Figure 10). The contribution of hydrogen bonding interactions of the new TCL analogues with the site residues as compared to the predominantly hydrophobic interaction of the reference compound suggests an elevated specificity of these virtual hits for the target (Figure 11).         a Ensemble average of the total energy E of the system (sum of potential energy E and kinetic energy E ); b Ensemble average of the potential energy E ; c IC was predicted by QSAR model (Table 3, Equation (2)) for the designed analogues. Experimental IC for the TCL11. * TCL11 is the co-crystallized ligand in the (PDB entry 1OOS) [21].  a Ensemble average of the total energy E of the system (sum of potential energy E and kinetic energy E ); b Ensemble average of the potential energy E ; c IC was predicted by QSAR model (Table 3, Equation (2)) for the designed analogues. Experimental IC for the TCL11. * TCL11 is the co-crystallized ligand in the (PDB entry 1OOS) [21].  a Ensemble average of the total energy E of the system (sum of potential energy E and kinetic energy E ); b Ensemble average of the potential energy E ; c IC was predicted by QSAR model (Table 3, Equation (2)) for the designed analogues. Experimental IC for the TCL11. * TCL11 is the co-crystallized ligand in the (PDB entry 1OOS) [21].
−114,415.0 −143,027.9 9.1 a Ensemble average of the total energy E tot of the system (sum of potential energy E pot and kinetic energy E kin ); b Ensemble average of the potential energy E pot ; c IC 50 was predicted by QSAR model (Table 3, Equation (2)

Binding Mode of TCL
Freundlich et al. reported the most comprehensive study of inhibition of PfENR by triclosan (TCL1) [21,23,27]. Characteristic intermolecular interactions of TCL1 and PfENR include -stacking interaction between the phenol moiety and the nicotinamide ring of NAD , hydrogen bonds between the hydroxyl group of TCL1 and -OH group of Tyr277 and the C 3 -OH group of the ribose ring of NAD . According to our analysis of the PfENR-TCLx complexes of predicted most potent inhibitors, these key interactions are conserved throughout this class of compounds and contribute also to predicted inhibitory potencies of novel triclosan derivatives. In addition, in the R -position of triclosan, a H-bond with the backbone of residue Ala217 and with the NAD cofactor was reported by Kumar et al. [31] (Figure 7). The most potent analogues TCL-60-08-08 and TCL-33-08-05 with the R -group ꞏꞏꞏCH 2 -cyclohexyl or the flexible ꞏꞏꞏoctyl do not have the capacity to form stacking interactions but can still well fill the hydrophobic pocket formed by residues Tyr267, Pro314, Ile369, Phe368, and Ala372. These R -groups, in combination with protonated amine or hydroxyl groups in the R -position and R -position that interact with the phosphate groups of NAD or with solvent, also contribute substantially to stabilisation of the PfENR-TCLx complexes. Thus, new TCL derivatives with potential to be developed into novel antimalarial agents that act on a new pharmacological target are proposed to medicinal chemistry laboratories for experimental verification.

Molecular Dynamics Simulations
Low deviations between the superimposition of active conformations of TCL inhibitors modelled by in situ modification of TCL11 and MM refinement and those obtained as ensemble averages over 500 snapshots from MD trajectories (Figure 13), suggest that the modelled complexes PfENR-TCLx are stable. Moreover, MD simulations confirmed the validity of the binding mode of new triclosan analogues generated from the crystal structure of PfENR-TCL11 [21].

Materials and Methods
The workflow describing consecutive steps of the combinatorial ligand design and the in silico process of novel TCL analogues against PfENR is presented in Scheme 1. Figure 13. Superimposition of the ligand active conformation from complexes refined by molecular mechanics (yellow carbons, red oxygens, blue nitrogens) and averaged active conformations resulting from MD simulations (purple carbon atoms).

Binding Mode of TCL
Freundlich et al. reported the most comprehensive study of inhibition of PfENR by triclosan (TCL1) [21,23,27]. Characteristic intermolecular interactions of TCL1 and Pf ENR include π-π stacking interaction between the phenol moiety and the nicotinamide ring of NAD + , hydrogen bonds between the hydroxyl group of TCL1 and -OH group of Tyr277 and the C 3 -OH group of the ribose ring of NAD + . According to our analysis of the Pf ENR-TCLx complexes of predicted most potent inhibitors, these key interactions are conserved throughout this class of compounds and contribute also to predicted inhibitory potencies of novel triclosan derivatives. In addition, in the R 3 -position of triclosan, a H-bond with the backbone of residue Ala217 and with the NAD + cofactor was reported by Kumar et al. [31] (Figure 7). The most potent analogues TCL-60-08-08 and TCL-33-08-05 with the R 1 -group ···CH 2 -cyclohexyl or the flexible ···octyl do not have the capacity to form stacking interactions but can still well fill the hydrophobic pocket formed by residues Tyr267, Pro314, Ile369, Phe368, and Ala372. These R 1 -groups, in combination with protonated amine or hydroxyl groups in the R 2 -position and R 3 -position that interact with the phosphate groups of NAD + or with solvent, also contribute substantially to stabilisation of the Pf ENR-TCLx complexes. Thus, new TCL derivatives with potential to be developed into novel antimalarial agents that act on a new pharmacological target are proposed to medicinal chemistry laboratories for experimental verification.

Molecular Dynamics Simulations
Low deviations between the superimposition of active conformations of TCL inhibitors modelled by in situ modification of TCL11 and MM refinement and those obtained as ensemble averages over 500 snapshots from MD trajectories (Figure 13), suggest that the modelled complexes Pf ENR-TCLx are stable. Moreover, MD simulations confirmed the validity of the binding mode of new triclosan analogues generated from the crystal structure of Pf ENR-TCL11 [21].

Materials and Methods
The workflow describing consecutive steps of the combinatorial ligand design and the in silico process of novel TCL analogues against Pf ENR is presented in Scheme 1.

Training and Validation Sets of TCL Inhibitors
The chemical structures and biological activities (IC exp 50 ) of training and validation sets of triclosan (TCL) based inhibitors of Pf ENR used in this study were taken from the literature [21,23,25,27]. The inhibitory potencies of these compounds cover a sufficiently broad range of half-maximal inhibitory concentrations (38 nM ≤ IC exp 50 ≤ 10 5 nM) to allow construction of a QSAR model. The training set (TS) contains 20 triclosan derivatives and the validation set (VS) contains 4 triclosan derivatives, taken from the refs. [21,23,25,27].

Model Building
Molecular models of the enzyme-inhibitor complexes Pf ENR-TCLx (E-I), free enzyme Pf ENR (E), and free inhibitors (I) were prepared from the high-resolution crystal structure of the reference complex Pf ENR-TCL11 (Protein Data Bank [43] entry code 2OOS at a resolution of 2.1 Å [21]) using Insight-II molecular modelling program [44]. The structures of the E and E-I complexes were considered at an appropriate physiologically relevant pH of 7 with capped neutral N-and C-terminal residues and all protonated and ionised residues charged. Crystallographic water molecules were included in the model. TS and VS were built into the structure of the reference complex Pf ENR-TCL11 by in situ replacement of the derivatised R groups of the TCL moiety. An exhaustive conformational search over all rotatable bonds of the replacing function group, coupled with a careful gradual energy optimisation of the modified inhibitor and the Pf ENR active site residues located in the vicinity of the inhibitor (within 5 Å distance), was employed to identify the low energy bound conformations of the modified inhibitor. The resulting low-energy structures of the E-I complexes were then carefully refined by minimisation of the whole protein. This procedure has previously been successfully used for model building of viral, bacterial, and protozoal protease inhibitor complex models and design of peptidomimetic, hydroxy naphthoic, and thymidine-based inhibitors [45][46][47][48][49][50].

Molecular Mechanics
The modelling of inhibitors, the Pf ENR and E-I complexes was carried out by molecular mechanics (MM) as previously described [47]. The TCL inhibitors, the Pf ENR enzyme, and their complex (E-I) were modelled in all-atom representation (parameters and charges) with the CFF force field. During all MM calculations, the dielectric constant was set to 4 to consider the dielectric shielding effect in proteins. The geometry optimisation process was performed gradually by relaxing the molecular structure through a sufficient number of steepest descent steps followed by iterative conjugate gradient cycles to reach the convergence criterion of 0.01 kcal·mol −1 ·Å −1 for the average gradient.

Conformational Search
The free inhibitor conformations were derived from their bound configurations in the E-I complexes by gradual relaxation to the nearest local energy minimum as previously described [46][47][48]. Monte Carlo search was carried out (for 50,000 iterations) over all rotatable bonds (excluding the rings) using Discovery Studio [34]. Two hundred unique stable conformations were generated for each TCL by random variation of the torsion angles of the last accepted conformer by ±15 • at 5000 K before energy minimisation with a dielectric constant ε = 80 to account for the dielectric screening effect of solvation. The lowest energy conformer was selected and minimised again at a dielectric constant of 4.

Solvation Gibbs Free Energies
The electrostatic component of Gibbs free energy (GFE) [50] that includes the effects of ionic strength by solving the nonlinear Poisson-Boltzmann equation [51,52] was computed by the DelPhi module in Discovery Studio [34] as previously described [46]. The implicit model of solvation represents the solvent by a continuous medium of high dielectric constant (ε o = 80) and the solute as charge distribution filling a low dielectric (ε i = 4) cavity with boundaries related to the molecular surface. The module numerically solves for the molecular electrostatic potential and reaction field around the solute using the finite difference method. DelPhi calculations were done on a (235 × 235 × 235) cubic lattice grid for the E-I complexes and free E and on a (65 × 65 × 65) grid for the free I. Full coulombic boundary conditions were used. Two subsequent focussing steps led to a similar final resolution of about 0.3 per grid unit at 70% filling of the grid by the solute. Physiological ionic strength of 0.145 mol·dm −3 , partial atomic charges and radii defined in the set of CFF force field parameters [34] and a probe sphere radius of 1.4 were used. The electrostatic component of the Poisson-Boltzmann solvation Gibbs free energy was calculated as the reaction field energy [53][54][55].

Calculation of the Entropic Contribution
The changes in vibrational entropy during the binding of the inhibitor to the enzyme (E) were calculated by the normal mode analysis of inhibitor vibrations using a simplified method based on Fischer et al. [56]. In this approach, vibrational analyses of the inhibitor bound at the active site of a 'frozen' E and of the low energy conformer of the free inhibitor are calculated for fully minimised structures.

Calculation of Binding Affinity and QSAR Model
The calculation of the binding affinity expressed as rGFE of complex formation has been fully described earlier [46][47][48].

Interaction Energy
The calculation of MM interaction energy (E int ) between enzyme residues and the inhibitor was performed as previously described [47,48].

Pharmacophore (PH4) Generation
According to the International Union of Pure and Applied Chemistry, a pharmacophore is "an ensemble of steric and electronic features that is necessary to ensure optimal supramolecular interactions with a specific biological target structure and to trigger (or block) its biological response" [57]. The evolution of the concept of pharmacophore in computer-assisted drug design has previously been previously reported [58][59][60]. The main task is to identify the highest spatial arrangement of common structural features shared by a group of molecules, making them recognisable by the active site of the biomolecular target and therefore characterising the affinity towards that target and consequently its biological activity. The first goal of the pharmacophore (PH4) is to estimate biological activity. Since the 3D Quantitative Structure Activity Relationships (3D-QSAR) Pharmacophore built here takes account of inhibitor-target interaction patterns, it is referred to as a structure-based model (differently from the model generated from a set of ligands without involving the target).
In this work, 3D-QSAR pharmacophore models were built using the Catalyst Hy-poGen algorithm [61] implemented in Discovery Studio [34] from the bound bioactive conformations of inhibitors taken from the QSAR models of E-I complexes. The E-I complexes were derived from the high resolution crystal structure of Pf ENR-TCL11 [21], so no new bioactive TCL conformations were generated during the process. Kurogi and Güner reported in a review the complete theory and successful applications of the PH4 generation process using the Catalyst software [62]. On the basis of the biological activities, HypoGen will optimise the pharmacophore hypotheses deduced from the training set inhibitors with the highest activity and omitting the least active molecules. The top scoring hypothesis of the pharmacophore will be built up in three main stages (constructive, subtractive and optimisation steps) from the set of most active inhibitors [62]. Inactive molecules serve for the definition of the excluded volume. The maximum number of five features allowed by the HypoGen algorithm was selected on the basis of the TCL scaffold and substituents during the pharmacophore generation, namely: hydrophobic aromatic (HYdAr), hydrophobic aliphatic (HYd), hydrogen-bond donor (HBD), hydrogen-bond acceptor (HBA) and ring aromatic (Ar). Many remaining features, among which positive charge (PC), positive ionisable (PIC) and negative charge (NC), negative ionisable (NIC) were not selected because of the neutral type of TCL molecules considered in this study. Additional features may be added if necessary [60]. Today, the concept of pharmacophore has been extended beyond the initial horizon of drug design to reach chemoinformatics and material-informatics QSPR as "Photovoltaphores" to identify metal-free dyes for dye-sensitized solar cells [63]. The adjustable parameters of the protocol were kept at their default values, except for the uncertainty of the activity, which was set to 1.25 instead of 3. Uncertainty is the range of the activity value with a default of <IC 50 /3, 3 IC 50 >. Our choice to bring the uncertainty interval for experimental activity to a relatively narrow <4 IC 50 /5, 5 IC 50 /4> takes into account the precision and homogeneity of the observed inhibitory activities measured by the same method in a single laboratory. Therefore, during the constructive step, the activity of the most active (MA) compound is multiplied by 1.25 instead of 3 to reach a value A = broadens the most active training set excessively. The subtractive step will consider as inactive those TCL compounds with activity 3.5 orders of magnitude beyond IC 50 (MA): IC 50 = 10 3.5 IC 50 (MA). Any pharmacophore that matches more than half of the set of inactive TCLs is removed. The default value of 3.5 is user adjustable if the range of experimental activities does not span 3.5 orders of magnitude [62]. For the optimisation phase, the hypothesis score will be improved for the remaining pharmacophore configurations by perturbations based on the estimated errors of the activity during the linear regression and complexity of each hypothesis. According to the law of parsimony (Occam's razor) "the simpler hypothesis that accurately estimates the activity is considered the best" [62].
During the generation of 10 pharmacophores, the number of missing features was set to 0 and the best pharmacophore was selected. The reliability of the generated pharmacophore models was assessed on the calculated cost parameters (in units of the number of bits required to generate them). The overall cost (or total cost) of a model is the sum of three components: (i) the weight, (ii) the error, and (iii) the configuration cost: Total cost = Cost(error) × CC(error) + Cost(weight) × CC(weight) + Cost(ConFigure) × CC(ConFigure) where CCs are cost coefficients in the Catalyst module of Discovery Studio [34]. Cost(error) and Cost(weight) are estimated through a Gaussian probability distribution: -Cost(error) increases as the Root Mean Square difference between the estimated and experimental activities for the training set compounds; -Cost(weight) increases in a Gaussian form as the feature weight deviates from an ideal value (2.0); - The configuration cost Cost(ConFigure) is a fixed cost depending on the complexity of the hypothesis space being optimised. It is equal to the entropy of the hypothesis space (log 2 P, P: the number of hypotheses initially created in a constructive phase that emerged through the subtractive one). In the standard HypoGen mode, its value should not exceed 17.
Details on costs and score calculations can be found in [49]. Two additional costs are calculated: the fixed cost and the null cost during the generation of the pharmacophore. They are theoretical costs helping to appreciate the significance of each PH4 hypothesis, since the closer the cost is to the fixed cost and the farther it is from the null cost the more significant it is statistically and the more likely to be a true correlation rather than a chance correlation.
Fixed cost represents the lowest possible cost of the perfect pharmacophore model that ideally fits the data (lower bound). Fixed costs are calculated by summing up the minimum possible error and weight cost plus the constant configuration cost. The null cost is another cost parameter; it stands for the maximum cost of a pharmacophore with no features and corresponds to the averaged activity data of the training set molecules (upper bound).

ADME Properties
Properties that determine the pharmacokinetic profile of a compound, in addition to the coefficient of octanol/water partitioning, aqueous solubility, blood/brain partition coefficient, Caco-2 cell permeability, serum protein binding, number of likely metabolic reactions and 18 more descriptors related to adsorption, distribution, metabolism and excretion (ADME properties) of the inhibitors, were computed by the QikProp programme [39] based on the methods of Jorgensen [40,41]. According to these methods, experimental results of more than 710 compounds were correlated with computed physicochemical descriptors, which resulted in an accurate prediction of the pharmacokinetic profile of considered new molecules. Drug likeness (#stars)-the number of property descriptors that fall outside the range of values determined for 95% of known drugs out of 24 selected descriptors computed by QikProp [39], was used as an additional ADME-related compound selection criterion.

Virtual Library Generation
The virtual library generation was performed as previously described [47]. Analogue model building was performed with the Molecular Operating Environment (MOE) programme [64]. The library of analogues was enumerated by attaching R groups (fragments, building blocks) to the MOE TCL scaffold using the Quasar CombiDesign module [64]. Chemical reagents considered in this study were taken from the directories of chemicals available from the commercial suppliers of chemicals [38]. Each analogue was built as a neutral molecule in the MOE program and its molecular geometry was refined by MM optimization using smart minimizer of Discovery Studio [34] with high convergence criteria (energy difference of 10 −4 kcal·mol −1 , R.M.S. displacement of 10 −5 Å) and a dielectric constant of 4 using class II consistent force field CFF [34] as described in Section 4.3.

ADME-Based Library Focussing
The drug-likeness selection criterion served to focus the initial virtual library as previously described [65].

Pharmacophore-Based Library Searching
The pharmacophore model (PH4) described in Section 4.8 and derived from the bound conformations of TCLs at the active site of Pf ENR served as a library search tool and in silico screening instrument as previously described [35]. The enumerated virtual library was further focused by using the ligand pharmacophore mapping protocol available from Discovery Studio [34]. Within this protocol, each generated conformer of the analogues was geometry optimised by means of the CFF force field for a maximum of 500 energy minimisation steps and subsequently aligned and mapped to the PH4 model to select the top-ranking overlaps. Twenty best-fitting conformers were saved and clustered into ten conformational families according to their mutual r.m.s.d by means of the Jarvis-Patrick complete linkage clustering method [66]. The best representative of each cluster was considered in the virtual screening of analogues.

Inhibitory Potency Prediction
The conformer with the best mapping on the PH4 pharmacophore in each cluster of the focused library subset was used for ∆∆G com calculation and IC pre 50 estimation by the complexation QSAR model as previously described [49]. The rGFE of formation of the E-I complex in water ∆∆G com was computed for each selected new analogue and then used to predict the Pf ENR inhibitory potencies (IC exp 50 ) of the focused virtual library of TCL analogues by inserting this parameter into the target-specific scoring function, Equation (2), Table 3. The scoring function, which is specific for the Pf ENR receptor: pIC exp 50 = a × ∆∆G com + b, was parameterised using the QSAR model described in Section 4.7.

Molecular Dynamics Simulations
The conformational stability of the co-crystallised ligand (TCL11) and the five selected top-ranking designed TCL analogues (TCL-33-08-05, TCL-58-01-01, TCL-58-16-08, TCL-59-01-01, TCL-60-01-01 obtained by in situ modification of the reference inhibitor TCL11 and refinement was evaluated by molecular dynamics (MD) simulations as described recently [67]. We have carried out 200 ns long simulations of these five explicitly solvated complexes in the NPT statistical ensemble (300 K, 1 bar) by using Desmond software [42]. A cubic periodic box with 10 Å buffer containing the E-I complex was filled with TIP3P water molecules and neutralised by adding the required number of counterions to reach electroneutrality (Table 9). During the simulation, OPLS4 force field [68], 1.5 fs integration step, and the coulombic interaction cut-off point of 14 Å, were used. The Nose-Hoover chains thermostat and Martyna-Tobias-Klein barostat methods were employed during the simulations [69,70]. After initial heating and relaxation, the productive simulation trajectory was recorded and analysed for ligand-receptor interactions every 400 ps, that is, 500 frames during each MD simulation. After MD calculations, the lowest total energy frame of each complex Pf ENR-inhibitor was minimised by molecular mechanics in 4 steps (as previously described [32]) using the OPLS4 force field in the Schrodinger software [42].

Conclusions
In this work, novel derivatives of TCL inhibiting ENR of Plasmodium falciparum that reach low nanomolar (one digit) inhibitory concentrations IC pre 50 have been designed (Table 8, Figure 12). A 3D-QSAR-based pharmacophore model of Pf ENR inhibition by a training set of 20 triclosan derivatives was used for in silico screening of a virtual combinatorial library of more than 33,000 triclosan analogues. More than 200 best virtual hits were modelled at the active site of Pf ENR using the crystal structure of Pf ENR-TCL11 complex [21] and their inhibitory potencies toward Pf ENR were predicted from computed rGFE of E-I complex formation, which correlated with the observed activities IC exp 50 of the training and validation set of TCLs. In silico screening, refinement of the structure of Pf ENR-TCLx complexes and prediction of IC pre 50 enabled identification of new TCL analogues with optimal substitutions in the R 1 , R 2 , and R 3 -positions of the TCL scaffold (Table 1). New triclosan analogues with predicted inhibitory potencies TCL-60-08-08 (IC pre 50 = 1.9 nM), TCL-33-08-05 (IC pre 50 = 4.1 nM), and TCL-33-16-08 with (IC pre 50 = 7.8 nM) display high predicted binding affinity to Pf ENR (Table 6), and favourable ADME-related properties (Table 7). Therefore, they could form new prospective candidates for reversible inhibitors aimed at the new validated pharmacological target of P. falciparum so far not connected with drug resistance to clinically used antimalarial therapeutics. Molecular dynamics simulations confirmed the stability of these enzyme-inhibitor complexes. Therefore, we can recommend these new molecules for synthesis and experimental testing of their inhibitory potency and antimalarial effect on the Pf ENR.