Synthesis and Hybrid SAR Property Modeling of Novel Cholinesterase Inhibitors †

A library of novel 4-{[(benzyloxy)carbonyl]amino}-2-hydroxybenzoic acid amides was designed and synthesized in order to provide potential acetyl- and butyrylcholinesterase (AChE/BChE) inhibitors; the in vitro inhibitory profile and selectivity index were specified. Benzyl(3-hydroxy-4-{[2-(trifluoromethoxy)phenyl]carbamoyl}phenyl)carbamate was the best AChE inhibitor with the inhibitory concentration of IC50 = 36.05 µM in the series, while benzyl{3-hydroxy-4-[(2-methoxyphenyl)carbamoyl]phenyl}-carbamate was the most potent BChE inhibitor (IC50 = 22.23 µM) with the highest selectivity for BChE (SI = 2.26). The cytotoxic effect was evaluated in vitro for promising AChE/BChE inhibitors. The newly synthesized adducts were subjected to the quantitative shape comparison with the generation of an averaged pharmacophore pattern. Noticeably, three pairs of fairly similar fluorine/bromine-containing compounds can potentially form the activity cliff that is manifested formally by high structure–activity landscape index (SALI) numerical values. The molecular docking study was conducted for the most potent AChE/BChE inhibitors, indicating that the hydrophobic interactions were overwhelmingly generated with Gln119, Asp70, Pro285, Thr120, and Trp82 aminoacid residues, while the hydrogen bond (HB)-donor ones were dominated with Thr120. π-stacking interactions were specified with the Trp82 aminoacid residue of chain A as well. Finally, the stability of chosen liganded enzymatic systems was assessed using the molecular dynamic simulations. An attempt was made to explain the noted differences of the selectivity index for the most potent molecules, especially those bearing unsubstituted and fluorinated methoxy group.


Introduction
As is known, cholinesterase (ChE) lyses choline-related esters that can serve as cholinergic neurotransmitters [1]. Roughly speaking, two groups of serine hydrolases are present in the human nervous system: acetyl-based (AChE) and butyryl-based (BChE) cholinesterases, respectively [2]. Unfortunately, the progressive degeneration of the central nervous system (CNS) with dysfunction of perception/cognition and dementia symptoms is characteristic for Alzheimer's disease (AD) [3]. Following the cholinergic hypothesis, the potential cholinesterase inhibitors (ChEIs) are postulated to reduce the psychosomatic AD disorders [4]. In fact, the competitive and non-competitive ChEIs are used as marketed drugs in the AD pharmacotherapy (e.g., galanthamine, rivastigmine, tacrine, donepezil) in order to diminish the disease symptoms [5][6][7][8].
Diversely substituted amide (-CONH-) or carbamate (-OCONH-) fragments can be incorporated as privileged structural motifs (scaffolds) in the design of the potential anti-ChE agents [9][10][11][12]. As a matter of fact, the amide/carbamate-based groups are present as building blocks or linkers in vast number of clinically implemented drugs [13][14][15]. Thus, the amide/carbamate-like anti-AD derivatives appear to be an attractive synthetic and modeling target as well.
The pharmacological property profiling and potency modeling strategies are commonly utilized to escort the synthetic efforts and nominate the prospective candidates in the hit→lead→seed→drug race [16]. SAR-guided (structure-activity) exploration of the descriptor-based feature space seems pivotal at the crucial decision-making phases of drug discovery [17]. Hence, in silico computation of the ADMET-tailored properties (Absorption, Distribution, Metabolism, Excretion, Toxicity) is expected to minimize the probability of drugs' late attrition [18]. As a matter of fact, a range of meaningful computer-supported approaches have been proposed to restrict the values of structural and/or physicochemical descriptors to the ADMET-friendly property space [19]. Basically, in silico-assisted approximation of the host-guest intermolecular recognition patterns can be dichotomized into the ligand-based and structure-related methodologies. In practice, the collaborative fusion of direct (receptor-dependent) and indirect (receptor-independent) procedures seems advisable in namely consensus approach for searching of promising drug molecule [20]. Conceptually, the mapping of the binding/active site in the pharmacophore-guided study is based on the simple tenet of substituent interchangeability and complementarity in the congeneric (structurally related) series of compounds. It is assumed that similar size, shape, and charge distributions should exert a comparable impact on the binding affinity and pharmacological profile, respectively [21]. Practically, the similarity-driven methodologies have been engaged in the comparable molecular shape analysis (e.g., CoMSA), where machine learning techniques are coupled with descriptor elimination/selection algorithms in order to approximate roughly the complex biological reality [22]. Despite some obstacles, the distance-related similarity investigations between structurally-alike molecules contribute significantly in the multidimensional quantitative (mD-Q) SAR practice [23]. Moreover, the numerical estimation of the ADMET-friendly features conjugated with desirable drug affinity (structure-activity landscape index, SALI) produces the SAR landscape, where the smooth/flat (homogenous) areas are striated with sharp (heterogeneous) activity cliffs [24]. Obviously, the clever management of SALI-related information can provide hints, which might be incorporated at the synthetic stage in order to modulate pharmacological effects and/or demolish unwanted side effects [25].
The relatedness of a ligand steric/electronic/lipophilic pattern can be also mediated by a complementary target surface using the protein-based docking procedures, where the (bio)effector-binding mode is determined according to the atomic coordinates of both host and guest molecules [26]. Despite drugs being basically engineered to fit tightly into the active interface of the enzyme/receptor, the binding goodness evaluation is continually questionable due to the lack of truly selective scoring functions [27]. Since the ligandenzyme interaction pattern is dependent on their conformations and relative orientations (poses), it seems advisable to investigate the guest-host system stability using molecular dynamic simulations (MDs) [28].
The primary objective of the presented study was the design, synthesis, and in vitro specification of the anti-AChE/BChE profile for the new set of ring-substituted benzyl [4-(arylcarbamoyl)phenyl-3-hydroxy]carbamates carrying different functionalized aryl side chains. Moreover, the structured in silico computation of ADMET-friendly property was performed to estimate the intermolecular similarity, lipophilicity, and inhibitory potential, respectively. The newly synthesized adducts were subjected to the quantitative shape comparison with the generation of an averaged pharmacophore pattern. The molecular docking approach was employed for the most potent AChE/BChE inhibitors in order to procure an additional comprehension of the binding mode. Finally, the stability of some fused ligand-enzyme systems was evaluated using the molecular dynamic simulations. An attempt was made to explain the observed variations of the selectivity index for the potent molecules.

Design and Synthesis
The synthesis of the carbamates of 4-aminosalicylanilides was carried out in two steps (see Scheme 1). The primary amino moiety of 4-aminosalicylic acid was protected by a reaction with benzyl chloroformate in an alkaline medium to form 4-{[(benzyloxy)carbonyl]amino}-2-hydroxybenzoic acid (CbzPAS) that was subsequently condensed with an appropriate substituted anilines using phosphorus trichloride in dry chlorobenzene under microwave conditions to give a series of investigated benzyl carbamates of 4-aminosalicylanilides. The lipophilicity of all the target compounds was investigated by means of RP-HPLC determination of capacity factors k with a subsequent calculation of log k. The retention times of individual compounds were determined under isocratic conditions with methanol as an organic modifier in the mobile phase using end-capped non-polar C18 stationary RP columns. All the results are shown in Figure 1 and Table 1.

Probability-Guided Pharmacophore Mapping
All synthesized ring-substituted benzyl [4-(arylcarbamoyl)phenyl-3-hydroxy]-carbamates 1-32 were tested to specify their in vitro inhibitory profile toward AChE and BChE targets and were compared with marketed drugs rivastigmine (RIV), galanthamine (GLT), and the pattern acid CbzPAS. Despite different mechanisms of action, RIV and GLT are the 2nd generation of ChEIs that pseudo-irreversibly inhibit AChE and BChE [8,14], respectively. Table 1 reports the IC 50 values (µM) that express the quantitative measure of inhibitor concentration needed to reduce the biological process by 50%.
Benzyl ( A comprehensive visualization tool for the investigation of SAR trends was provided in the form of the structure-activity landscape indexes (SALI), where the planar image of molecular similarity is augmented with molecular activity profile [29]. Basically, the detection of the continuity areas and/or activity cliffs depends largely on the accessibility of structurally related compounds (chemotypes) with noticeable variations in the activity profile. The distribution of the Tanimoto coefficients was analyzed for the triangular T 32×32 matrix revealing that the majority of compounds are characterized by pretty high similarity metric (0.85 ≥ T ≥ 0.75). Interestingly, the potent compounds 2 and 15 are marked as fairly similar with T = 0.93. Obviously, for structurally related molecules (e.g., stereoisomers T→1), SALI→ infinity; therefore it is recommended to replace such values by the largest SALI values [30]. The symmetrical grayscaled heatmap of SALI values for the congeneric series of analyzed compounds is illustrated in Figure 2a, where axes correspond to molecules sorted according to the increasing pBChE inhibitory potential (∆pBChE = 1.2), respectively. Accordingly, the white spot of the heatmap represent the highest numerical values of SALI parameters, while the black ones specify the minimal ones. In fact, molecule 15 is characterized by low values of SALI indexes with the remaining compounds in the dataset, while the most potent compound 2 is accompanied by gray pots of the inactive molecules 3, 4 (positional isomers) and 1 (unsubstituted analogue), respectively. The structural modifications that considerably affect the molecular activity profile are illustrated as the neighborhood plot in Figure 2b, where the structural relatedness between pairs of compounds is shown in the function of the activity difference and colorcoded according to SALI values. Noticeably, three pairs of fairy similar fluorine/brominecontaining compounds (17 vs. 19 and 28, 26 vs. 32) can potentially form the activity cliff that is manifested formally by high SALI numerical values (see Figure 2b). In consequence, further dense samplings of the indicated SAR-variations (T > 0.95 & ∆pBChE > 1) seem reasonable for the investigated set of BChE inhibitors. The crucial distributions of the electronic/steric/lipophilic properties of the guest-host composition can be specified by the systematic sampling of the functional group relatedness with the introduction of the pharmacophore concept in QSAR studies [31]. In reality, the spatial arrangement of pharmacophoric features is determined for the homogeneous ensemble of molecules that share the common structural pattern (chemotype), because the enhancement of the SAR applicability domain (AD) for a non-congeneric set of structures is a rather elusive and enigmatic operation. Moreover, the prediction of the modeled property is firmly dependent on the molecule training and test separation. On the other hand, there are no strictly defined rules for splitting the set of compounds into training/test subfamilies in order to maintain the model robustness as well as its predictive power. Hence, the recurrent and interchangeable training/test subset division was proposed for the probability-driven pharmacophore mapping based on the stochastic model validation (SMV) algorithm [32]. Despite the CPU-intense (Central Processing Unit) SAR calculations, the whole pool of systematically produced training/test subpopulations (ratio 3:1) was investigated in CoMSA modeling of the BChE inhibition potency (C 8 32 ≈ 10 7 models). The frequency distribution in the test set for preferable statistic parameters of model robustness (q 2 cv > 0.6 and q 2 test > 0.5) indicted that the potent compound 2 is noticeably overrepresented in contrast to molecule 15. Surprisingly, the preferential selection of the active compounds (2 and 6) and inactive ones (3, 13, 23 and 26) that are basically orthoand meta-substituted analogues resulted in the generation of the robust models with the acceptable predictive power for the test set. Finally, the selection-oriented pharmacophore pattern was generated based on the robust models following the principles of IVE (iterative variable elimination-partial least squares)-PLS algorithm that is described elsewhere [33]. The graphical representation of the descriptor-based areas that contribute (un)favorably into CoMSA models is shown in Figure 3. The relative contribution of the surface/charge descriptors is weighted by the corresponding regression coefficient indicating the regions of the positive (bright color) and negative (dark color) impact on the inhibitory potency (see Figure 3a). Moreover, the four possible combinations of the charge and the mean regression coefficients are introduced in Figure 3b.
Noticeably, a diverse set of steric and electrostatic features was indicated on the averaged surface for the enzyme-selective ligands; therefore, the direct translation of the pharmacophore-based areas of the space into the complementary pseudoreceptor model that potentially harbors putative inhibitors is a fairy complicated task. The spatial pattern plotted in Figure 3a indicates the favorable steric contributions of areas that spread uniformly over the 4-{[(benzyloxy)carbonyl]amino}-2-hydroxybenzoic acid skeleton and aniline-like substituent as well. Not surprisingly, the positively charged areas of hydrogen atoms directly attached to nitrogen atoms and the negatively charged oxygen atoms of carbonyl groups were depicted as the favorable contributors to the inhibitory potency of the investigated compounds (see Figure 3b). The specified atoms can potentially form the hydrogen bonds with the hypothetical receptor aminoacid residue in the enzyme active site. Unfortunately, the provided pharmacophoric pattern does not explain the observed variations in selectivity of molecule 2 and 15 toward A/BChE enzymes. In fact, no pharmacophoric areas were specifically pointed out as dominant ones to increase the inhibitory potency; however, the consensus 3D-QSAR modeling provides the spatial map of steric/electronic features of the ligand-site composition.

Similarity-Driven Property Assesment
The distance-oriented property evaluation was conducted using the Principal Component Analysis (PCA) and Hierarchical Clustering Analysis (HCA) on the pool of 2819 descriptors derived from Dragon 6.0 software. The descriptor-based data were organized into a matrix X 32×2819 with rows presenting molecules (objects) and columns representing numerical variables (parameters), respectively. The resulting dataset was centered and standardized, since the calculated parameters varied considerably. The number of significant principal components (PCs) was determined according to the percentage of the modeled variance. The first four PCs accounted for 77.6% of the total data variance, while the first three PCs described 72.42%. The respective scoreplots with the projection of molecules 1-32 and color-coded with the corresponding molecular weight (MW) on plane PC1 vs. PC2 are presented in Figure 4a   . Interestingly, the lower weighted molecules (MW < 400) are located in the negative area of the plane (PC1 < 0 and PC2 < 0).
Moreover, the exploratory analysis of the descriptor-driven data structure was performed by tracing the (dis)similarities between objects (molecules) in the multidimensional variable space in order to reveal its clustering tendency. Due to its hierachical nature, the HCA method basically leads to sub-optimal clustering of objects that is mainly dependent on the procedure employed for clusters' linkage [34]. The interpretability of the extracted data structure is not obvious in the high-dimensional variable space; therefore, the findings are presented in the form of dendogram constructed on the Euclidean-based distance and the Ward linkage algorithm, where the OX axis presents the sequence of objects/parameters and the OY axis specifies the dissimilarity [35]. The dendogram shown in Figure 5a proved our previous PCA findings (see Figure 4a), where at least three sub-clusters (A,B,C) can be distinguished-molecules 12-15 and 27 (blue lines) differ noticeably from the rest of the objects (red lines), respectively. The HCA dendogram was augmented with a colorcoded map of the experimental activity and lipophilicity data in the logarithmic scale as presented in Figure 5b. As a matter of fact, a simultaneous interpretation of the dendogram objects (sorted according to the Ward linkage method) and the colored map of the empirical data revealed that there is no clear relationship between the investigated data. Roughly speaking, the molecules grouped in cluster D are characterized by relatively high values of the log k parameter when compared to the remaining objects. From the user's perspective, the non-trivial question arises of how to single out the best-suited method or combination of methods for clogP estimation of new compounds. In order to investigate the lipophilic profile of the analyzed molecules, the experimentally specified log k parameters were compared with the calculated lipophilicity coefficients (clogP) estimated using a set of in silico predictors, e.g., clogP, Molinspiration, Osiris, HyperChem 7.0, Sybyl-X, MarvinSketch 15, ACD/ChemSketch 2015, Dragon 6.0, Kowwin, XlogP3, ChemDraw, and ACD/Percepta programs (see Table S1 in the Supplementary Materials). Moreover, the selected physicochemical properties of rivastigmine, galanthamine, and compound 2 were predicted by ACD/Percepta 2012 (see Table S2 in Supplementary Materials). The corresponding clogP estimators deduced by the group of alternative programs were cross-compared with the experimental potency values and inter-correlated with each other, as illustrated in Figure 6. Relatively good correlations (r = 0.77) between the generated clogP values using Sybyl-X program and the empirical log k were recorded. The averaged clogP values over the set of programs produced the value of r = 0.73, whereas the median resulted in r = 0.70 with the experimental data. On the whole, a pretty high inter-correlation was observed among the programs for lipophlicity specification (see Figure 6); however, some variations in clogP values resulted probably from different computational algorithms (atom/fragment-or descriptor-based) implemented in the software and/or training data engaged at the training stage-models are as good as the modeling data used. Additionally, the iterative variable elimination procedure (IVE-PLS) was applied on the integrated clogP matrix (X 32 × 13 ) and log k parameter indicating Sybyl-X and XlogP3 estimators as significant contributors to the linear quantitative structure-property relationship QSPR model in namely consensus lipophilicity determination [36]. Consequently, the averaged values of the indicated clogP estimators were correlated with log k values resulting in r = 0.80, because not only the best inter-correlated logP estimators were chosen. Similarly, the same programs were specified with the stepwise procedure implemented in Matlab environment.

Docking and Molecular Dynamics Simulations
The pharmacophore-driven (ligand-based) SAR approaches can be coupled with the site-directed (structure-based) docking procedures that position a potential bio-effector molecule within a protein binding/active. The complementary host-guest binding mode is deduced based on the target spatial arrangement of atoms with the engagement of feature/descriptor-matching algorithms, where the ligand property space is correspondingly mapped to the receptor steric, electrostatic, and/or lipophilic patterns. Unfortunately, it is still not clear how to directly relate the pharmacological or ADMET effects to enthalpically and/or enthropically favorable ligand-receptor modes and scoring descriptors. On the other hand, in silico attempts to reconstruct the drug-protein interactions using the molecular docking are commonly accepted as a complimentary approach to the traditional mD-QSAR ligand-based methods [37].
The structural geometries of human acetyl-/butyrylcholinesterase enzymes (A/BChE) in a complex with the commercially accessible drug molecules (e.g., galanthamine and rivastigmine) are being currently intensively scrutinized. The crystallographic structure of BChE determined at a sophisticated resolution of 2.6 Å in the ligand-containing state (holo) with the rivastigmine (RIV) analogue [3-[(1~{R})-1-(dimethylamino)ethyl]-4oxidanyl-phenyl]~{N}-ethyl-~{N}-methyl-carbamate was retrieved from the Protein Data Bank repository (PDB code: 6eyf). Subsequently, the marketed drugs and the potential anti-BChE population (IC 50 < 150) were docked in the active site of the enzymatic chain A using AutoDock Vina software in order to collate the binding mode of the analyzed analogues with the drug-enzyme interacting pattern [38]. The resulting molecular conformations and orientations (poses) of the potentially active compounds revealed mainly three types of non-binding interactions: the hydrogen bond (HB), hydrophobic, and π-stacking interactions, respectively. The graphical illustrations of the relevant ligand-enzyme interactions were produced for the entire population of the active compounds (IC 50 < 150) using Schrödinger Maestro software and Protein-Ligand Interaction Profiler (PLIP) and compared with the commercial drug molecules (GLT and RIV) [39]. As a matter of fact, the hydrophobic interactions were overwhelmingly generated (70%) with Gln119 (20%), Asp70 (13.7%), Pro285 (13.7%), Thr120 (12.7%), and Trp82 (9.5%) aminoacid residues, while the HB-donor ones were dominated (50%) with Thr120, which confirms our previous findings [1,2]. π-stacking interactions were specified with Trp82 aminoacid residue of chain A as well. The planar (2D) and spatial (3D) binding pattern distribution for marketed drugs and the most active anti-BChE molecule 2 are presented in Figures 7 and 8, respectively. The hydroxyl group of Thr120 seems to be the key structural motif in the formation of hydrogen bond (HBD) with hydroxyl oxygen of molecule 2 (see Figure 8c). The π-stacking interaction was revealed between the substituted anilide-based ring and Trp82, respectively. Unfortunately, the replacement of hydrogen atoms in the methoxy group of molecule 2 with fluorine in the trifluoromethoxy substituent in compound 15 does not change noticeably the enzyme-ligand binding pattern; however the additional hydrogen bond (HBA) is generated using the oxygen (acceptor) of Asn68 carbonyl group and nitrogen (donor) of molecule 15. The molecular electrostatic potentials (MESP) as well as the lipophilic potentials determined on Connolly surface of molecule 2 and 15 were determined. In fact, the interrelation was reported previously between electrostatic and lipophilic potentials, because the molecular charge distribution affects the lipophilic property [40]. Noticeably, the electron-rich substituents correspond well to higher lipophilic potential values, as illustrated in Figure 9. Not surprisingly, the variations of charge distribution in molecule 2 and 15 have a visible impact on the lipophilic patterns, especially in the regions of -OCH 3 and -OCF 3 groups as indicated by brow colors (see Figure 9).  The molecular dynamic simulations (MDs) can supply fruitful knowledge regarding the host-guest complex stability as well as the potential mechanism of ligand action in the active site of the enzyme [41]. The liganded state (holo) of the BChE enzyme docked with the most potent molecule 2 was applied as the input structure in the MD simulations. Taking into account the technical feasibility (computational time and resource constraints), the CPUintense MD calculations were conducted with the involvement of the flexible ligand and the rigid target according to the procedure implemented in Sybyl-X 2.2.1 program. Contrarily, the internal macromolecular motions (atom trajectories) and the system energy variations were investigated based on the reduced molecular dynamics (RedMD) algorithm. Accordingly, in the RedMD strategy, each aminoacid/nucleotide is expressed by a singular spherical particle ('bead') centered on a Cα or P atom with the mass that corresponds to the total mass of a given aminoacid or nucleotide [42][43][44]. The Berendsen algorithm was engaged in RedMDs on the NVT ensemble (constant number, volume, temperature) in order to control temperature. The initial atom velocities were assigned at 293K; however, the system was then heated, and the corresponding trajectories were specified at the temperature of 300K. The total length of the RedMDs was arranged to 20 ns with a 10 fs time step, respectively. The enzyme fluctuations of the total energy during the RedMD simulations are illustrated in Figure 10a. In fact, relatively small variations of the root mean square deviation (RMSD) values indicated the enzyme stability, as shown in Figure 10b. Noticeably, the side chains of compound 2 and 15 have an important impact on the resulting variations of RMSD values recorded for heavy atoms during MDs, as depicted in Figure 11. It seems that trifluoromethoxy substituent of molecule 15 is characterized by lower spatial flexibility compared to the methoxy group of compound 2. Obviously, -OCF 3 is considerably heavier and can generate greater steric hindrance in comparison to -OCH 3 . The most electronegative element fluorine can interact with the surrounding aminoacid residues that stiffen the ligand-enzyme system, as mirrored in Figure 11.

In Vitro Cell Viability Assay
To evaluate the potential cytotoxic effect, the preliminary in vitro screening of the antiproliferative activity of the most effective AChE/BChE inhibiting compounds 2 and 28 was performed. The human monocytic leukemia THP-1 cell line was selected for this purpose, as we described previously [45]. Nor 2 neither 28 show negative effect on relative cell viability up to the concentration of 30 µM. The relative cell viability was around 100% at this concentration. It makes compounds 2 and 28 interesting for further pharmacological evaluation and allows them to become lead compounds.  [46]: 4-Aminosylicylic acid (6.0 g, 39.2 mM) was disolved in methanole (64 mL) at ambient temperature, and benzyl chloroformate (6.8 mL, 47.6 mM) was added dropwise. The reaction mixture was stirred at room temperature for the time period of 24 h. Then, the solvent was removed to dryness under reduced pressure. Crude product was dissolved in the mixture of ethyl acetate (200 mL) and 1M hydrochloric acid (200 mL). Layers were divided. The water layer was additionally extracted with ethyl acetate (3 × 50 mL). Afterwards, organic layers were collected and dried with magnesium sulfate. Then, the solvent was removed to dryness under reduced pressure. Yield 71%; Mp 219-220 • C; IR (ATR, cm −1 ) : 1727, 1623, 1589, 1538, 1440, 1383, 1307, 1270,  1249, 1220, 1204, 1185, 1165, 1104, 1056, 1028 4-{[(Benzyloxy)carbonyl]amino}-2-hydroxybenzoic acid (0.6 mmol; 0.2 g) was suspended in dry chlorobenzene (6 mL) at ambient temperature and phosphorus trichloride (0.3 mM, 0.5 eq.), and the corresponding substituted aniline (0.6 mM, 1 eq.) was added dropwise. The reaction mixture was transferred to the microwave reactor, where the synthesis was performed (30 min, 130 • C). Then, the mixture was cooled to 40 • C, and then the solvent was removed to dryness under reduced pressure. The residue was washed with hydrochloride acid and water. The crude product was recrystallized from ethanol.

Determination of Lipophilicity by HPLC
The HPLC separation system Agilent 1200 series equipped with a DAD SL (Agilent Technologies) was used. A chromatographic column Symmetry ® C18 5 µm, 4.6 × 250 mm, Part No. W21751W016 (Waters Corp., Milford, MA, USA) was applied. The HPLC separation process was monitored by the ChemStation for LC 3D chromatography software (Agilent Technologies). Isocratic elution by a mixture of MeOH p.a. (72%) and H 2 O-HPLC Mili-Q grade (28%) as a mobile phase was used for the determination of the capacity factor k. The total flow of the column was 1.0 mL/min, injection 20 µL, column temperature 40 • C, and sample temperature 10 • C. The detection wavelength of 210 nm was chosen. A KI methanolic solution was used for the determination of the dead times (t D ). Retention times (t R ) were measured in minutes. The capacity factors k were calculated according to the formula k = (t R − t D )/t D , where t R is the retention time of the solute and t D is the dead time obtained using an unretained analyte. Each experiment was repeated three times. The log k values of individual compounds are shown in Table 1.

Evaluating In Vitro AChE and BChE-Inhibition Potencies
The ability of all the prepared compounds to inhibit AChE from electric eel (Electrophorus electricus) and BChE from equine serum (both purchased from Sigma, St. Louis, MO, USA) was determined in vitro using a modified Ellman's method, as described in detail elsewhere [47]. This is a simple, rapid, and direct method to determine the SH and -S-Sgroup content in proteins. The enzyme activity is measured indirectly by quantifying the concentration of the 5-thio-2-nitrobenzoic acid ion formed in the reaction between the thiol reagent 5,5'-dithiobis-2nitrobenzoic acid (DTNB) and thiocholine, which is a product of acetylthiocholine hydrolysis catalyzed by cholinesterases [48][49][50]. The effectiveness of the tested inhibitors is expressed as the IC 50 values.
The enzyme activity in the final reaction mixture (2 mL) was 0.2 U/mL, the concentration of the substrate (acetylthiocholine or butyrythiocholine) was 40 µM, and the concentration of DTNB was 0.1 mM for all reactions. All reactions were carried out at 25 • C in the presence of phosphate-buffered saline (0.1 M, pH 7.4) in a glass cuvette with a 1 cm optical path. All the studied inhibitors were dissolved in DMSO (concentration 0.01 M) and then diluted in demineralized water as necessary. For all tested compounds and standards, rivastigmine and galanthamine in at least four different concentrations of inhibitor in the final reaction mixture were used. All measurements were carried out in triplicate, and the average values of the reaction rate (v 0 -reaction in the absence of the inhibitor, v i -reaction in the presence of the inhibitor) were used for construction of the dependence v 0 /v i vs. concentration of the inhibitor. The IC 50 values were calculated from an obtained equation of the regression curve for y = 2 (resulting from the definition of IC 50 ) [14,15,25]. The results are presented in Table 1.

In Vitro Viability Assay
Human monocytic leukemia THP-1 cells were used for in vitro antiproliferative assays. Cells were obtained from the European Collection of Cell Cultures (ECACC, Salisbury, UK) and routinely cultured in Roswell Park Memorial Institute RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS), 2% L-glutamine, 1% penicillin and streptomycin (all from Sigma-Aldrich) at 37 • C with 5% CO 2 . Cells were passaged twice a week. The effect of the compounds on cell viability was determined using a Cell Counting Kit-8 (CCK-8, 2-(2-methoxy-4-nitrophenyl)-3-(4-nitrophenyl)-5-(2,4-disulfophenyl)-2H-tetrazolium, monosodium salt, Sigma) according to the manufacturer's instructions. The compounds were dissolved in DMSO and added to the cell suspension in the RPMI 1640 medium without FBS (5 × 10 5 cells/mL). The maximum concentration of DMSO (Sigma) in the assays never exceeded 0.1% (v/v). Subsequently, the cells were incubated at 37 • C with 5% CO 2 for 24 h. After this incubation period, the relative cell viability (the ratio between cells treated with compounds and cells treated only with DMSO) was determined by CCK-8 kit.

Model Building and Molecular Modeling
CACTVS/csed and CORINA editors were applied to generate molecule structural model with its 3D geometry. The data format conversion was performed using an Open-Babel (inter)change file format converter. The Sybyl-X 2.0/Certara package installed on a DELL workstation with Ubuntu 20.10 operating system was engaged in order to conduct the molecular modeling simulations. The MAXMIN2 module implemented in Sybyl-X was used to initially optimize the compound spatial geometry with the standard Tripos force field (POWELL conjugate gradient algorithm) with a 0.01 kcal/mol energy gradient convergence criterion. The electrostatic potential values were calculated using Gasteiger-Hückel method. One 18-ordered atom trial alignment was produced on the most active compound 2 (active analogue approach) with an FIT procedure to cover the entire bonding topology in the maximal common structure (MCS). SONNIA software was employed to simulate 10 × 10 to 30 × 30 self-organizing maps (SOMs) with a winning distance in the range from 0.2 to 2.0 in CoMSA analysis. Spatial coordinates of the molecular surfaces and the corresponding potential values were used as input to Kohonen SOM network to generate a 2D map of the electrostatic potential (MEP) for the set of superimposed molecules. The produced maps were reshaped into a 100-to 900-element vector subjected to the PLS method implemented in the MATLAB environment.

Principal Component and Hierarchical Custering Analysis
The human-friendly 2D/3D illustrations of the compound distribution in the experimentalbased (FCS) and virtual-derived (VCS) molecular space might be displayed by Principal Component Analysis (PCA). PCA is a linear projection method that can be engaged to model multidimensional data with a relatively small number of so-called principal components (scores and loadings) produced to maximize the description of variance within the input data.
The PCA model with f principal components for a data matrix X can be calculated as follows in Equation (1): where X is a data matrix with m objects and n variables, T is the score matrix with dimensions (m × f ), P T is a transposed matrix of loadings with dimensions (f × n), and E is a matrix of the residual variance (m × n) not explained by the first f principal components.
On the whole, the first few principal components (PCs) frequently describe sufficiently data variance and reveal the groups of objects. Hierarchical Clustering Analysis (HCA) enables us to examine the (dis)similarities between objects in the variable space. Hence, the similarity measure as well as the way of resulting sub-clusters linkage should be specified a priori. The produced findings are illustrated as a dendogram, where the OX axis presents the indices of the clustered objects, while the OY one corresponds to the linkage distances between two objects linked, respectively. Moreover, the visualization method can be applied to the empirical data sorted according to the order of objects with the generation of the color-coded variable maps. A mutual interpretation of objects sorted with the Ward linkage method and the color-coded variable maps facilitate the (dis)similarity assessment of objects in terms of the input parameters.

Theoretical Lipophilicity Evaluation
A number of freely/commercially available in silico estimators can be employed to calculate the theoretical partition coefficients (clogP) as follows: AlogPS-procedure proposed by Tetko et al. based on atom-type electrotopologicalstate (E-state) indices and neural networks (NN); milogP-method invented by Molinspiration for practical logP calculations of almost all organic molecules as a sum of fragment-based contributions and correction factors; ClogP-fragment-based procedure for predicting lipophilicity based on structuredependent correction values retrieved from Hansch and Leo's database that is implemented in Sybyl/Centara software; HyperChem logP-an atom-additive method that approaches lipophilicity using the individual atomic contribution based on procedure proposed by Ghose, Prichett, and Crippen; MarvinSketch logP-the overall lipophilicity of a molecule is composed of the contributing values of its atom types that were redefined to accommodate electron delocalization and contributions of ionic forms; ChemSketch logP-a comprehensive fragment-based algorithm with high-quality models retrieved using the empirical data. Well-characterized logP contributions have been compiled for atoms, structural fragments, and intramolecular interactions provided for more than 12 × 10 3 experimental logP values; Dragon AlogP-the statistical estimators of the Ghose-Crippen-Viswanadhan model were calculated on the basis of known experimental logP for the training set of 8364 compounds. The overall estimation of the lipophilic atomic-based constant is evaluated with the contribution of 115 atom types; Dragon MlogP-the theoretical partition coefficient includes VdW volume and Moriguchi polar parameters as correction factors. A regression MlogP model is based on 13 structural parameters, which were evaluated on the training group of 1230 organic molecules; Kowwin-estimates the log octanol-water partition coefficient of chemicals using the atom/fragment contribution approach; XlogP3-an atom-additive procedure with well-defined correction factors that employs an optimized atom typing approach calibrated on a big training set; OSIRIS clogP-in-house approach based on the cumulative sum of atom contributions calculated for more than 5000 compounds with experimentally determined logP values as training set. Predicting engine distinguishes 369 atom types; ChemBio clogP-the algorithms for prediction of partition coefficient based on a training set of compound provides coverage for a broad chemical space; Percepta clogP-based on >12 × 10 3 of experimental logP values with the algorithm that uses the principal of isolating carbons.
The redundant descriptors of QSAR/QSPR investigations were selected and eradicated using the modified version the uninformative variable elimination (UVE-PLS) method, namely iterative variable elimination (IVE-PLS). Briefly, the entire algorithm is composed of the following stages: (1) standard PLS analysis with leave-one-out crossvalidation LOO-CV to evaluate the performance of the PLS model; (2) elimination of the matrix column with the lowest abs(mean(b)/std(b)) value; (3) standard PLS analysis of the new matrix without the column eliminated in (2); (3) iterative repetition of (1)-(3) to maximize q 2 cv value.

Similarity-Based Activity Landscape Index
The numerical profiling of similarity-related structure-activity landscape index (SALI) can be quantitatively expressed according to the following Equation (2) where A x and A y are the activity profiles for the x-th and y-th molecule and sim(x,y) is the pair-wise similarity measure. The Tanimoto coefficient was engaged for the fingerprintbased similarity estimation, where the structural pair-wise molecular relatedness is calculated as follows; in Equation (3): T(x, y) = n xy n x + n y − n xy (3) where n xy is the number of bits set into 1 shared in the fingerprint of the molecule x and y, n x is the number of bits set into 1 in the molecule x, n y is the number of bits set into 1 in the molecule y, respectively.

Ligand/Structure-Based Analysis
The selection-driven surface analysis of a structural space can be performed using the set of descriptors that grab information about the spatial pharmacophoric pattern. In practice, some machine learning methods conjugated with weighting and selecting procedures were engaged to indicate the minimal ensemble of pharmacophoric features that are potentially important in description of guest-host interactions. The comparative molecular surface analysis (CoMSA) was applied in order to meaningfully compare the shape and charges variations generated on the molecular surface of the receptor and ligand, respectively. Briefly speaking, a single layer of neurons is arranged in a 2D plane with welldefined topology to produce self-organized maps (SOMs). The spatially closed objects are located in the proximal neurons of the square map in the process of SOMs adaptation to the input data. In consequence, a 2D image of the property space is produced, where structurally related molecules are placed in neighboring neurons. The electrostatic/steric features that are potentially important in the ligand-receptor recognition and complementarity phenomena can be specified using the iterative variable selection approaches. The backward column eradication was recurrently repeated until the optimal number of variables included within the model was accomplished-the moment that the q 2 cv deterioration specified the set of potentially relevant columns. The cumulative sum of the common columns for all of the investigated BChE models was calculated and normalized to a range of (0,1). The crystallographic structure of butyrylcholinesterase determined using X-ray diffraction at 2.6 Å resolution co-crystalized with a rivastigmine analogue (liganded state) was downloaded from the PDB repository (PDB code: 6eyf). Apart from the RIV-based analogue ([3-[(1~{R})-1-(dimethylamino)ethyl]-4-oxidanyl-phenyl]~{N}-ethyl-~{N}-methylcarbamate), all remaining heteroatoms (including water molecules) were eradicated prior to docking in the AutoDock Vina programe. Initially, the ligand/enzyme structures were prepared in the pdbqt file format with the calculated Gasteiger charges. The grid box was centered on the central atom of the RIV analogue (BY2). In AutoDock Vina, docking simulations in different poses (default nine) were generated progressively from a single conformer (an energy-minimized molecule). Then, the resulting molecular conformations and orientations with the preferred torsion angles and the rotatable bonds were evaluated by the unitedatom (UA) scoring function. VMD, PyMol, Schrödinger Maestro graphical viewers, and Protein-Ligand Interaction Profiler (PLIP) were employed to illustrate the foreseen 2D/3D binding modes, respectively.
Sybyl-X 2.0/Certara modeling software was used to conduct the dynamic simulations of the ligand-enzyme systems. The host-guest structures were analyzed and fixed by the Prepare Protein Structure module. The Berendsen algorithm was used in MD simulations on the NVT ensemble in order to control temperature. The initial atom velocities were assigned at 293K; however, the system was then heated and the corresponding trajectories were determined at the temperature of 300K. The total length of MDs was arranged to 20 ns with 10 fs time step, respectively. The RedMDStream 1.0 program with coarsegrained molecular dynamics models was engaged for the macromolecular simulations in the same conditions.

Conclusions
In summary:

3.
The distribution of the Tanimoto coefficients was analyzed for the triangular T 32×32 matrix. Interestingly, the potent molecules 2 and 15 are marked as fairly similar (T = 0.93). In fact, molecule 15 is characterized by low values of SALI indexes with the remaining compounds in the dataset, while the most potent molecule 2 is accompanied by gray pots of the inactive molecules 3, 4 (positional isomers), and 1 (unsubstituted analogue), respectively. Noticeably, three pairs of fairy similar fluorine/bromine-containing compounds (17 vs. 19 and 28, 26 vs. 32) can potentially form the activity cliff that is manifested formally by high SALI numerical values.

4.
The preferential selection of the active molecules (2 and 6) and inactive ones (3, 13, 23, and 26) that are basically orthoand meta-substituted analogues resulted in the generation of the robust models with the acceptable predictive power for the test set. 5.
The molecular docking approach was engaged for the most potent AChE/BChE inhibitors in order to get comprehensive knowledge of the binding mode. The hy-drophobic interactions were overwhelmingly generated with Gln119, Asp70, Pro285, Thr120, and Trp82 aminoacid residues, while the HB-donor ones were dominated with Thr120. π-stacking interactions were specified with Trp82 aminoacid residue of chain A as well. 6. The stability of some fused ligand-enzyme systems were evaluated using the molecular dynamic simulations. The trifluoromethoxy substituent of molecule 15 is characterized by lower spatial flexibility compared to methoxy group of compound 2. Obviously, -OCF 3 is considerably heavier and can generate greater steric hindrance in comparison to -OCH 3 , as the most electronegative element fluorine can interact with the surrounding aminoacid residues, which stiffens the ligand-enzyme system.