Multi-combined QSAR, molecular docking, molecular dynamics simulation, and ADMET of Flavonoid derivatives as potent cholinesterase inhibitors

Abstract In searching for a new and efficient therapeutic agent against Alzheimer’s disease, a Quantitative structure-activity relationship (QSAR) was derived for 45 Flavonoid derivatives recently synthesized and evaluated as cholinesterase inhibitors. The multiple linear regression method (MLR) was adopted to develop an adequate mathematical model that describes the relationship between a variety of molecular descriptors of the studied compounds and their biological activities (cholinesterase inhibitors). Golbraikh and Tropsha criteria were applied to verify the validity of the built model. The built MLR model was statistically reliable, robust, and predictive (R2 = 0.801, Q2cv = 0.876, R2test = 0.824). Dreiding energy and Molar Refractivity were the major factors that govern the Anti-cholinesterase activity. These results were further exploited to design a new series of Flavonoid derivatives with higher Anti-cholinesterase activities than the existing ones. Thereafter, molecular docking and molecular dynamic studies were performed to predict the binding types of the designed compounds and to investigate their stability at the active site of the Butyrylcholinestérase BuChE protein. The negative and low binding affinity calculated for all designed compounds shows that designed compound 1 has a favorable affinity for the 4TPK. Moreover, molecular dynamics simulation studies confirmed the stability of designed compound 1 in the active pocket of 4TPK over 100 ns. Finally, the ADMET analysis was incorporated to analyze the pharmacokinetics and toxicity parameters. The designed compounds were found to meet the ADMET descriptor criteria at an acceptable level having respectable intestinal permeability and water solubility and can reach the intended destinations. Communicated by Ramaswamy H. Sarma


Introduction
Neurodegenerative diseases are progressive dysfunctions of the nervous system caused by a deterioration of the functioning of nerve cells.Alzheimer's disease is among the complex neurodegenerative diseases with many cognitive and neuropsychiatric manifestations (Mehta et al., 2012).This disease is characterized by a disruption of brain tissue that leads to the progressive loss of cognitive functions especially memory, language, visuospatial performance, and executive functions (Varikasuvu et al., 2019;Wan et al., 2020).Important behavioral symptoms include restlessness, apathy, irritability, anxiety, illusions, hallucinations, and disinhibition (Crapser et al., 2020).Therefore, the manifestations of this disease lead to progressive disability and eventual incapacity.
The alteration of cholinergic function is of crucial importance in Alzheimer's disease, especially in brain regions that deal with learning, memory, behavior and emotional responses, including neocortex and hippocampus (Takousis et al., 2019).Cerebral atrophy is the most obvious clinical discovery in Alzheimer's disease where levels of acetylcholine (ACh) (Anand & Singh, 2013), a neurotransmitter responsible for conducting electrical impulses from one nerve cell to another nerve cell, are decreased due to its rapid hydrolysis by acetylcholinesterase enzyme (AChE) (Amin et al., 2021;Vecchio et al., 2021).This enzyme produces non-cholinergic secondary functions that include promoting beta-amyloid deposition in the form of senile plaques/neurofibrillary tangles in the brain of affected individuals (Ahmed et al., 2021;Marucci et al., 2021).The initiation and progression of Alzheimer's disease can be due to beta-amyloid deposition.
Butyrylcholinesterase (BuChE) is an enzyme that regulates cholinergic neurotransmission using the hydrolysis of Ach (Zhou & Huang, 2022).The development of Alzheimer's disease shows an increase in BuChE in the most affected areas of the brain such as the hippocampus and the temporal cortex (Kumar et al., 2020).The increase of BuChE allows to beta-amyloid aggregation in the early stages of senile plaque formation.Though, BuChE is mainly located in peripheral tissues, including plasma and a very small amount is present in the brain region (Jasiecki & Wasa ˛g, 2019;Taha et al., 2021).
Therefore, AChE and BuChE are considered critical targets for the successful management of Alzheimer's disease by increasing the availability of ACh in brain regions and decreasing beta-amyloid deposition (Kashyap et al., 2020).However, BuChE is still associated with the neuropathological structures of Alzheimer's disease, neuritic plaques, and neurofibrillary entanglement.The decrease in AChE with the progression of Alzheimer's disease can produce reducing yields to increase acetylcholine levels in the brain which decrease since to cholinergic neurotoxicity.In contrast, the selective inhibition of BuChE increases the levels of acetylcholine in the brain.This prompted the suggestion that BuChE may be a better target than AChE for the treatment and development of new drugs for Alzheimer's disease (Okello & Mather, 2020;Stellenboom, 2019;Uddin et al., 2021).
Grace to the development of computer science and the growth of intensive parallel computing, in particular, molecular modeling has become a real challenge.The research and synthesis of new chemical and biochemical compounds are now often associated with a molecular modeling study.Various theoretical chemistry methods were investigated to determine the physical or chemical properties of molecules, including quantum chemistry methods for accurately determining the electronic properties of molecules and molecular mechanics methods which are based on empirical parameters to determine the structural and physicochemical parameters (Bannwarth et al., 2021;Verma & Truhlar, 2020).The determination of the physicochemical parameters was carried out using the QSAR study which consists of estimates of the effects of a variation in molecular structure on biological activity (Abramenko et al., 2020;Adhikari et al., 2020).Numerous studies show the effectiveness of the QSAR method in the estimation of biological molecules and enzymes for chronic diseases (Ghanei-Nasab et al., 2021;Guo et al., 2019;Hammoudan et al., 2021;Kurniawan et al., 2019;Oyinloye et al., 2021;Ponzoni et al., 2019).
In recent years, the pharmaceutical industry was turned to molecular modeling as a new research method to predict the properties and activities of drug molecules before they are synthesized.In this context, the present work focuses on the application of QSAR molecular modeling methods to predict expected biological activities in new molecules for Alzheimer's disease.The aim was to develop and evaluate the QSAR model for the modeling of a new series of flavonoid derivatives with higher anti-cholinesterase activities.The molecular dynamic model was used to predict the anti-cholinesterase activity of newly designed flavonoid compounds.A molecular anchoring study was performed to determine the binding mode responsible for the stability of the molecules at the active site of the BuChE protein.

Data set
A QSAR study was conducted on 45 flavonoid compounds that had previously been evaluated for their anti-cholinesterase properties (Yingying & Hongwei, 2019).The experimental activity values, IC50, representing the concentration required to inhibit 50% replication of Alzheimer's disease for the studied compounds, were converted to the negative logarithm of IC50 (pIC50 ¼ -log (IC50)).The chemical compositions of the examined compounds, along with their experimental activity levels (pIC50), are listed in Table S1.Furthermore, to investigate the features that govern anti-cholinesterase activity, several molecular descriptors were calculated using three modeling programs ChemOffice (ChemOffice, n.d.), ChemSketch (ACD/ChemSketch Freeware, 2016), and MarvinSketch (MarvinSketch, n.d.).

Statistical methods
The molecular descriptors generated for the studied compounds were analyzed using principal component analysis (PCA) (Jolliffe & Cadima, 2016) to reduce their dimensionality without losing the initial information.The resulting dataset was then randomly divided into training and test sets in a 5:1 ratio.Furthermore, the multiple linear regression (MLR) techniques in the XLSTAT software (XLSTAT version 2019, n.d.) were employed to develop a linear model that describes the relationship between the calculated descriptors and anti-cholinesterase activity.This software is user-friendly and easily accessible.

Model validation
The statistical parameters proposed by Golbraikh and Tropsha (Golbraikh & Tropsha, 2002;Tropsha, 2010) including the coefficient of determination (R 2 ), the adjusted coefficient for degrees of freedom ðR 2 adj ), mean square error of the model (MSE), Fischer coefficient (F), variance inflation factor (VIF), cross-validation coefficient ðQ 2 cv ) and coefficient of determination for the test set ðR 2 test Þ were calculated and analyzed to assess the fit, stability, and predictive power of the resulting QSAR model.Furthermore, the Y-randomization test was performed to verify the robustness of the obtained model.
According to the Organization for Economic Cooperation and Development (OECD) Principle 3 (37), a QSAR model is valid only within its applicability domain (AD).Therefore, it is essential to determine the built model's AD before making new predictions.AD is a technique for identifying compounds outside the bounds of the constructed QSAR model, and it recognizes outliers in the compounds of the training set.The QSAR models' applicability domain (AD) can be defined in several ways (Tunkel et al., 2005), but the most typical method is to calculate the leverage values, where xi is the row vector of the query compound descriptor, and X is the n � k À 1 ð Þ matrix of k model and the descriptor values for n compounds in the training set.The superscript t denotes matrix/vector transposition (Qin et al., 2017).In this study, we used the Williams plot, and the leverage threshold in this plot is h et al., 2005), where n is the number of compounds in the training set, and k is the number of model descriptors.The domain of applicability is defined as falling within a square area within three standard deviations (in this study, x ¼ ±3; 'rules of three sigmas') from the mean value.

Molecular docking
Molecular docking is a well-established method used to investigate the binding patterns of a ligand and receptor, including the amino acids involved in biological activity (Atatreh et al., 2019;Madhavi Sastry et al., 2013).Before proceeding with molecular docking, we downloaded the structure of BuChE from the RCSB database (Oyinloye et al., 2021) using the PDB ID 4TPK.The 4TPK entry corresponds to the complex of human butyrylcholinesterase with a co-crystallized ligand (Ref).The protein was prepared using Discovery Studio software (Dassault Syst� emes, n.d.), by removing all water molecules and the co-crystallized ligand, and then adding polar hydrogens, Kollman charges, and missing atoms.The active site of BuChE was defined by the sphere containing the co-crystallized ligand.We employed AutoDock Tools (AutoDock, n.d.) and Vina (AutoDock Vina, 2022) for molecular docking, using the AUTOGRID algorithm to define the three-dimensional grid and determine the binding energy between the ligand and the receptor (Morris et al., 1998).The grid size was set to x ¼ 22 Å, y ¼ 32 Å, and z ¼ 30 Å, with a distance of 0.375 Å between grid points.The center of the grid was set to the active site of the receptor, which had coordinates of (x ¼ 4.02 Å, y ¼ 10.76 Å, and z ¼ 13.86 Å).Finally, the docking results obtained using AutoDock Tools and Vina were visualized using Discovery Studio software.

Molecular dynamics simulation
The protein-ligand complex was performed using molecular dynamics using the GROMACS 2019.3 (B� o et al., 2020;Naz et al.).The GROMOS96 54a7 force field was used for protein, and the topology for ligands was generated using the PRODRG server (Verma et al., 2016).All the complexes were solvated using simple point charge (SPC) water molecules in a rectangular box (Wu et al., 2006), and the systems' net charges were neutralized by adding 3 Cl ions.The steepest descent method achieved both minimum energy and maximum force, with Fmax set at 1000 kJ/mol/nm.Two successive 1000 ps simulations using canonical NVT and isobaric NPT ensembles were used to equilibrate the system at 300 Kelvin and a pressure of 1 bar (Paci & Marchi, 1996;Tabti et al., 2023).Following this, 100 ns of molecular dynamics simulations were run for each molecule.The output trajectories were produced, and the files were evaluated to grasp the protein behavior better.The graphs were created with the Ubuntu version of the gnu plot plotting tool.

Binding energy calculation by MM-PBSA
To determine the binding free energy of 4TPK-Ref and 4TPKdesigned molecule 1 complexes and validate the ligandprotein interaction investigation, the researchers employed the molecular mechanic/Poisson-Boltzmann surface area (MM-PBSA) methods.The binding energy was calculated by the 'g_mmpbsa' packages using GROMACS, and it includes potential energy and polar and non-polar solvation energies.The binding energy calculations in this method were calculated using the following equation: binding represents the total binding energy of the complex, while the G receptor and G ligand represent the binding energy of the free receptor and unbounded ligand, respectively.

Molecular descriptors
As previously mentioned, a QSAR study was conducted on a series of 45 Flavonoid derivatives to establish a quantitative relationship between their structural descriptors and the studied biological activity.To investigate the primary structural descriptors that govern anti-cholinesterase activity, we calculated 37 descriptors for each molecule (Table 1) using ChemSketch (ACD/ChemSketch Freeware, 2016), ChemOffice (ChemOffice, n.d.), and MarvinSketch software (MarvinSketch, n.d.).The values of these 37 descriptors are presented in Tables S2-S4.Furthermore, we employed the principal component analysis (PCA) method to identify correlations between the computed descriptors and eliminate those with high correlations.As a result, the number of descriptors was reduced to twenty, which are listed in Table S5.

Data splitting and model development
Eighty percent of the molecules in the dataset were chosen randomly as the training set, while the remaining twenty percent were used as the test set.Then, the multiple linear regression (MLR) methods was applied to build a mathematical model that explains the anti-cholinesterase activity variance.Several MLR models were generated, although, models that did not meet the Golbraikh and Tropsha criteria (Golbraikh & Tropsha, 2002;Tropsha, 2010) or the OECD principles (OECD, 2014) were categorically eliminated.Thus, the top fifty MLR models in terms of R 2 coefficients of determination, indicating that the selected descriptors explain the majority of Anti-cholinesterase variance, the highest predictive ability proved by their high R 2 test values, highest cross-validation coefficient Q 2 cv and the lowest values of MSE (Tables S6-S8).The R 2 adj values for all the created models are relatively close to R 2 , suggesting that the number of descriptors used in the models is appropriate and overfitting is not a concern (B� o et al., 2020).The low MSE values further support this conclusion.Additionally, the high Q 2 cv values in cross-validation indicate that the models are statistically robust.The strong predictive power of the models is evident from the high R 2 test values.Moreover, the developed models meet the cutoff points for several statistical metrics proposed by Golbraikh, Tropsha, and other experts.Table S7 presents the statistical parameters for the selected QSAR model Equations, all of which have high correlation coefficients (R 2 > 0.6; R 2 adj > 0.6), high leave-one-out cross-validation determination coefficients (Q 2 CV > 0.6), high external validation determination coefficients (R 2 test > 0.6), and low mean squared error.Hence, these models demonstrate high internal and external predictive power.The high correlation (R 2 ) between the predicted and observed activities of the investigated compounds suggests the reliability of the models.Furthermore, the high Q 2 CV and R 2 test values confirm the robustness of the 50 models and their ability to predict the anti-BuChE activity of new compounds within their applicability domain.Notably, all the built models have R 2 values that are close to R 2 adj values, indicating that they are not overfitted.Model 28 emerges as the best MLR model based on several statistical metrics ( R 2 and Figure 1 shows the predicted (pIC50 (Pred)) and observed (pIC50 (Obs)) activities of model 28, respectively.Assuming the null hypothesis (i.e.no effect of explanatory factors), we would expect a risk of less than 0.01% because the p-value is less than 0.0001.Therefore, we can confidently state that the chosen molecular descriptors provide a significant amount of information.Moreover, the selected descriptors exhibit no multicollinearity, and the resulting model is stable as indicated by the VIF values of the tree descriptors, which are less than 10 (Verma et al., 2016) (VIF ¼ 6.441 and 6.442 for DE and MR, respectively).

Y-randomization test
To confirm the robustness of the model, we repeated all computations with the randomized activities of the training set compounds (Y randomization test).We generated 100 randomized trials using stand-alone QSAR tools available online at Cheminformatics Tools and DTC Lab Tools.As shown in Table 3, none of the randomized trials could replicate the original model.The average values of R Rand , R 2 Rand and Q 2 cv ðRandÞ are 0.226, 0.065, and À 0.128, respectively, the value of c?? 2 is 0.775 and all the new random models have significantly low R 2 Rand and Q 2 cv ðRandÞ values, These results confirm that the developed QSAR model is robust and was not built by chance.

Applicability domain (AD)
To evaluate the applicability domain (AD) of model 28, we conducted a leverage analysis and plotted the normalized

Design of new compounds
To design new compounds with higher anti-cholinesterase activity, the primary descriptors governing the activity were explored using model 28.It was found that MR (Molar refractivity) positively influenced activity, while DE (Dreiding energy) had a negative impact.To evaluate the importance of each descriptor on the variance of pIC50, standardized coefficients or t-test values were necessary.A higher absolute t-test value indicates a more influential descriptor.In our model, the t-test values for MR and DE were 8.167 and À 4.300, respectively, suggesting that MR has a greater influence on the activity than DE.By analyzing the descriptors in the selected QSAR model, it is possible to determine the factors connected to the anti-cholinesterase activity.Therefore, the selected descriptors can be interpreted as follows: -Dreiding energy (DE) is related to the 3D structure of the molecule, and it is computed using the Dreiding force field by adding the energy components.This force field energy can be used to differentiate between several conformations of the same molecule based on their stability.In this work, the dreiding energy was calculated using MarvinSketch software (Table S5).In general, the force fields are mainly applied to compute the interaction energy between the protein and the ligand as pair-wise interaction potentials consisting of van der Waals and electrostatic interactions, in addition to H-bond energy between the ligand and the enzyme (Baassi et al., 2023;Psimadas et al., 2012;Soufi et al., 2023;Tuenter et al., 2017).In Eq. ( 1), DE is characterized by a negative sign, indicating that there is an inverse relationship between anti-cholinesterase activity and this descriptor.In the literature several studies were investigates the Dreiding energy (G� omez-Jeria et al., 2021;Jackson et al., 2008;Jain & Jadhav, 2013).Therefore, Neha Kansal was investigate the QSAR 3D study performed using molecular field analysis to identify some of the important regions of the inhibitor that have electrostatic interactions and specific sterile with three different kinases including KDR, cKit and FLT3 using the least genetic partial square (G/PLS) method (Kansal et al., 2008).These AMF models for the three receptors revealed that ligand molecules should have a diaryl fraction urea with the 3-Amino pyrazolo pyridine or 3-amino indazole nucleus.However, this 3D-QSAR study examines the essential structural features, which can be used to synthesize new and powerful multitarget RTK inhibitors for cancer prevention.Though, Tabish Equbal was studied using the 3D-QSAR model of the inhibitory activity of FTase has been developed from sterile and electrostatic descriptors to study the substitution requirements for the receptor-favorable drug and to predict the relative inhibitory activities of 37 analogues of SCH 66336 (Equbal et al., 2008).A stable and statistically significant model with a high correlation coefficient was achieved.A significant predictive capability of the observed model for the molecules of the external test set supports that the derived model can be used for the design of new inhibitors.This 3D-QSAR study examines the necessary structural features, which can be used for modifications in analogues SCH 66336 to achieve better FTase inhibitory activity.
-Molecular refractivity (MR) is an important property in several chemical and pharmacological processes.It is a measure of the substance's total polarizability per mole and is dependent on the molecule's lipophilicity and polarizability, providing information about the molecule's volume and cohesion.The Molar refractivity can be determined as a function of temperature, refractive index, and pressure.In Eq. ( 1), MR is characterized by a positive sign, indicating that there is a direct proportionality between MR and anti-cholinesterase activity.
Therefore, the design of new molecules with a high anticholinesterase activity requires: -A high MR value by adding stronger withdrawing electron groups.-A low DE value by adding stronger electron donor groups.
Considering the above description, 23 new molecules with improved activity values than the studied ones (Table 4) were designed, and their activities pIC50 were predicted using model 28 (Eq.1).
The results indicate that all the compounds that were designed have high predicted pIC50 values and fall within the applicability domain of model 28, as indicated by the h � value being less than 0.208.These findings demonstrate the successful collaboration between theoretical and pharmacological approaches.

Molecular docking
All of the newly created compounds were subjected to molecular docking studies within the binding site of the human butyrylcholinesterase (BuChE) protein, using the crystal structure with RCSB Code: 4TPK as a reference (Brus et al., 2014).The results, as summarized in Table 5, revealed that all of the designed compounds had docking score values that ranged from À 9.6 to À 11.8 kcal/mol.The co-crystallized ligand formed two Carbon Hydrogen bonds interactions with Gly A: 116, and Glu A: 197, one Pi-Anion bond with Asp A: 70, it exhibited also hydrophobic interactions with Ile A: 69, Gly A: 116, Trp A: 231, Phe A: 329, Leu A: 286, and Tyr A: 332 residues (Tables 6, 7).According to Table 5, ligand 1 exhibited the highest binding affinity towards the target protein (BuChE).To gain insights into its mode of interaction with the protein, the involved interactions were analyzed.As shown in Table 7, ligand 1 exhibited a complex set of molecular interactions that allowed it to form a stable complex with the target protein.Indeed, it turned out that it implicates one Conventional Hydrogen Bond with His A: 438, as well as a Halogen bond interaction with Pro A: 285, and a Pi-Pi T-shaped interaction with Trp A: 231 and Phe A: 329, Pi-Alkyl interactions with Trp A: 82, Tyr A: 332, Phe A: 329, His A: 438, Ala A: 328, and Leu A: 286.Furthermore, an Amide-Pi Stacked interaction was observed with Gly A: 116, along with a Pi-Pi Stacked interaction with Tyr A: 332.The presence of hydrophobic and conventional hydrogen bond interactions in the stable complex indicates that ligand 1 is a more potent inhibitor of BuChE.Finaly,we selected the best docked complex (ligand 1-BuChE) with higher binding affinity and the best interaction for molecular dynamics to deep understanding their stability during 100 ns.

Validation docking protocol
The docking protocol has been validated by the 're-docking' technique.Co-crystal structure BuChE of PDB-4TPK was chosen and was redocked back into their active site.The redocking complex was then superimposed on the BuChE co-crystallized ligand using PyMOL.The RMSD value of the selected target was found to be 1.50 A˚ Figure 3, this low RMSD value indicates that the molecular docking protocol could be reliable.Therefore, the ligands could interact with the same amino acid residues as those reported in the cocrystallized BuChE (Tabti et al., 2022).

Molecular dynamics simulation
Molecular dynamics modeling was used to assess the effect of the most active experimental (Ref) and designed compound 1 on the structure of Butyrylcholinesterase (BuChE) And their stability in the enzyme's binding pocket (Lu et al., 2019).Accordingly, Gromacs2019.3 was used to conduct a molecular dynamics simulation (100 ns) of the complexes 4TPK-Ref and 4TPK-designed molecule 1 (Figure 4).For that, different characteristics such as root mean square deviation (RMSD), Root-mean-square fluctuation (RMSF), and radius of gyration (Rg) were examined on the trajectories' corresponding data files.The stability of complex protein-ligand structures is typically evaluated using RMSD, RMSF, and Rg plots.

Conformational protein
The RMSD was computed during the simulation to establish the overall stability of the selected systems and considered the calculated value as the primary criterion for measuring the system's convergence (Baammi et al., 2022(Baammi et al., , 2023)).The RMSD values of the 4TPK -designed molecule 1 complexes are less (0.20 nm) than the 4TPK -Ref complex (0.27 nm) over the entire simulations confirming the higher stability of 4TPK-designed molecule 1.Indeed, the backbone RMSDs plot analysis showed increasing trends with increasing RMSD values between 0 and 70 in the case of 4TPK -Ref and 4TPK -designed molecule 1, suggesting that the Ref and the compound 1 were becoming accustomed to a new conformation inside the binding pocket, after that, the plateau continued and finally settled at 0.30 ± 0.02 nm and 0.18 ± 0.02 nm respectively which does not exceed the 0.3 nm threshold (Figure 4A).

Residue mobility analysis
The root mean square fluctuation (RMSF) method was used to investigate the effect of the Ref and designed molecule 1 binding on the dynamics of the target protein during 100 ns binding on the flexible structure of the protein as well as the behavior of essential amino acids (Khan et al., 2021).
Table 7. 3D and 2D binding interactions of the 2 complex compounds.

Radius of gyration analysis
We examined how the compactness of the protein's structure is being modified as a consequence of the linkage to various ligands (En-Nahli et al., 2022;Stockwell & Thornton, 2006).For this goal, we calculated the radius of gyration (Rg) as a function of time.Any ligand with a high enough Rg value will be more likely to be flexible, making it unstable.On the other hand, lower Rg values point to confirm that it is dense and closely packed confirmation (Aljuaid et al., 2022).
The average Rg values of 4TPK -Ref and 4TPK-compound 1 were 2.29 ± 0.05 nm and 2.28 ± 0.00, respectively.Moreover, based on the graph (Figure 4C), it appears that the Rg of 4TPK -designed compound 1 had tighter packing than 4TPK-Ref.

Hydrogen bonds
The hydrogen bond is an essential feature that determines binding affinity and contributes to the binding relationship between ligands and proteins (Schiebel et al., 2018).In drug discovery, it is also in charge of drug specificity, metabolization, and adsorption.To confirm the stability of all complex compounds, the hydrogen bonds, and pairs within 0.35 nm between protein-Ref, protein-designed molecule 1 were estimated in a solvent environment during the MD simulations (Figure 4D, E).A comparison of the Ref, the molecule 1 made more hydrogen bonds with the target protein (4TPK), indicating a strong binding activity.

MMPBSA analysis
The MM-PBSA approach was used to calculate the binding free energy (DG) between 4TPK-Ref and 4TPK-designed molecule 1 complex for the whole trajectories which are displayed in Table 8.The interaction between a ligand and protein is shown in the form of binding energy.The final binding energy is a cumulative sum of van der Wall, electrostatic, polar solvation, and SASA energies.The average binding free energy of the complexes including Ref, and designed molecule 1 was À 56.599, and À 45.809 kJ/mol respectively.As shown in Table 8, the newly designed compound 1 had the lowest binding energy than the FDA-approved drug indicating that molecule 1 reduce the binding strength.Except for the polar solvation energy, all other forms of energy contributed favorably to the interaction between different molecules and 4TPK which supports docking results.

Bioavailability and ADMET parameters screening
Many medication development phase failures are brought on by unfavorable pharmacokinetic and toxicity issues.The drug discovery process would greatly benefit us if these problems could be resolved at an early stage.To investigate the pharmacokinetic and toxicity parameters for 22 drugs, the ADMET (absorption, distribution, metabolism, elimination, and toxicity) properties were utilized.These entities were tested in humans for their ADMET properties, aqueous solubility, blood-brain barrier (BBB), CYP2D6 binding, intestinal absorption, hepatotoxicity, and plasma protein binding (PPB).The pharmacokinetic profiles of all the anti-cholinesterase of a series of Flavonoid derivatives compounds under investigation were predicted in Discovery Studio 2.5 using ADMET descriptors (Accelrys, San Diego, CA, USA).A collection of rules/keys that determine threshold ADMET features for the chemical structure of the molecules based on the given drug information are used in the module's six mathematical models to forecast properties quantitatively (Paci & Marchi, 1996).Table 9 presents the analysis findings together with a biplot (Figure 5).Human intestinal absorption (HIA) following oral delivery is predicted by ADMET absorption.The 95% and 99% confidence ellipses in the ADMET PSA 2D and ADMET AlogP98 plane, respectively, define the absorption levels of the HIA and BBB penetration model (Pires et al., 2015).The areas between the ellipses are those that are anticipated to contain well-absorbed chemicals.The compound's hydrophilicity can be estimated from the AlogP98 value.PSA is another crucial factor that affects medication bioavailability because drugs with PSA ˂ 140 Å 2 can be passively absorbed, leading to high oral bioavailability (Lever et al., 2016;Roy et al., 2012).The 1,2,4,5,6,11,12,13,21,and 4, 5, 6, 7, 8, 9, 11, 12, and 18 had AlogP98 values � 5 (Table 9).These physicochemical factors are connected to tolerable intestinal permeability and water solubility, which are the initial stages of oral bioavailability.Most drugs have unknown BBB penetration levels, as shown in (Table 9).The prospective medications' water solubility is a key factor in their bioavailability, with compounds 4, 5, and 6 having an acceptable level (level 2) while the other compounds have extremely low levels.Absorption of all molecules is moderate to low.According to the CYP2D6 inhibition predictions, these 22 hits are not CYP2D6 inhibitors (Table 9).This shows that all substances are efficiently digested during Phase-I metabolism and that no substance causes harmful medication interactions.ADMET hepatotoxicity forecasts possible human hepatotoxicity for a wide spectrum of structurally varied chemicals.Drugs' bioavailability may be decreased if they bind to plasma proteins (Du et al., 2016;Lu et al., 2019).All 22 hits had weak PPB characteristics, according to the ADMET plasma protein binding property prediction (binding rates ˂90%).
In conclusion, Compounds 1, 2, 4, 5, 11, and 13 were discovered to meet the ADMET descriptor requirements at an acceptable level with appropriate water solubility and intestinal permeability based on the in silico ADMET study (Table 9) (Figure 5) and can accomplish the required targets.

Conclusions
This work was concerned with modeling the quantitative structure-activity relationship for the development of QSAR models capable of predicting the anti-cholinesterase activity of a series of 45 flavonoid derivatives.To check the validity of the built model, two criteria Tropsha and Golbraikh were applied.A high correlation coefficient (R 2 ¼ 0.801, Q2cv ¼ 0.876, R 2 test ¼ 0.824) was detected confirming the reliability, robustness, and predictive of the multiple linear regression model.The anti-cholinesterase activity was influenced by molar refractivity and binding energy.The negative binding energy values of all designed compounds (-9.9 to À 11.6 Kcal/mol) show a favorable affinity and stability between the ligands and the protein Ligand-4TPK complex.Molecular dynamic simulation detects the stability of compound 13 in the active pocket of the 4TPK over 100 ns.
ADMET analysis was incorporated to analyze pharmacokinetics and toxicity parameters.The criteria of the ADMET descriptors were adapted to the designed compounds at an acceptable level, respecting water solubility and intestinal permeability, which could reach the intended destinations.
According to the results, the QSAR simulation method shows a very high efficiency in the estimation of Alzheimer's disease drugs which reduces the time and cost of syntheses of these drugs which can be applied in the pharmaceutical industries.

Figure 1 .
Figure 1.(A) Correlations of observed and predicted activities calculated by the (model 28).(B) Contribution of model descriptors (model 28).
of ±3 standard deviations (x ¼ 3).Furthermore, none of the training or test set compounds are outliers as all are within the AD.Only one molecule (molecule 24) in the training set has a leverage value (hi) above the limit (h � ¼ 0.208), while molecule 1 in the training set is the only compound predicted to fall outside the standard deviation limit (± 3r).

Figure 2 .
Figure 2. Williams plot of the residual normalized against the leverage for the MLR model (with h � ¼ 0.208 and residual Limits ¼ ± 3); form the samples in black color and test the samples in red color).

Figure
Figure 4B illustrates the RMSFs for the 4TPK -Ref and 4TPKdesigned molecule 1 complexes.Compared to 4TPK -Ref, the complex related to designed compound 1 exhibits weaker fluctuations, not exceeding 0.3 nm.Furthermore, the average RMSF values for 4TPK -Ref and 4TPK -molecule1 were 0.12 ± 0.05 nm and 0.07 ± 0.04 nm, respectively.These results indicate that the binding of compound 1 to the 4TPK caused No important variations in protein conformation.

Figure 4 .
Figure 4.The results of the molecular dynamics study: (A) the time evolution of the target protein's backbone, (B) The comparative RMSF values for the target protein with the Ref and designed molecule 1, (C) The comparative radius of gyration values for the target protein with the Ref and designed molecule 1, (D, E) the comparative Hydrogen bonds and pairs within 0.35 nm values for the target protein with the Ref and designed molecule 1.
22 chemicals have 99% confidence levels for human intestinal absorption, as illustrated in Figure 5.The remaining compounds fall outside the filter of the ADME model ellipses, indicating their poor intestinal absorption and BBB penetration ability.The model states that a chemical must meet the requirements (PSA < 140A˚2 and AlogP98 � 5) to have optimal cell permeability (Falc� on-Cano et al., 2020; Roy et al., 2012).All compounds had polar surface areas that met the PSA < 140A˚2 criterion.Considering the AlogP98 criteria, compounds 1, 2,

Figure 5 .
Figure 5. Prediction of drug absorption for 22 anti-cholinesterase activity of a series of Flavonoid derivatives hits.Discovery Studio 2.1 (Accelrys.San Diego.CA) ADMET descriptors.2D polar surface area (PSA 2D) in � A˚2 for each compound is plotted against their corresponding calculated atom-type partition coefficient (ALogP98).

Table 1 .
List of computed descriptors for the studied molecules.

Table 2 .
Observed and predicted activities, residuals, and leverage for model 28.
residuals and the leverage threshold values, where h � ¼0.208(h � ¼3 � (k þ 1)/n), k ¼ 3,and n ¼ 36.Any predicted pIC50 value outside the compounds entering the AD on which the model was built cannot be considered reliable.If the leverage value (h) of a compound is greater than the threshold (h � ), it suggests that the compound has a significant influence on the model.The outliers have a hi > h � and normalized residuals higher than 3. Based on the results shown in Figure2, all compounds in the training and test sets fall within the range

Table 3 .
Results of the Y-randomization test for model 28.

Table 4 .
Structures, predicted activity, and leverage values of newly designed compounds.

Table 6 .
The 2D interactions between the protein BuChE and the ligand 1 and co-crystallized ligand (Ref).

Table 5 .
The binding affinity of designed compounds (kcal/mol).

Table 9 .
In silico ADMET properties of the 22 hits.

Table 8 .
MM-PBSA calculations of binding free energy for Reference drug and designed molecule 1.