Design of novel anti-cancer drugs targeting TRKs inhibitors based 3D QSAR, molecular docking and molecular dynamics simulation

Abstract Tropomyosin receptor kinase (TRK) enzymes are responsible for different types of tumors caused by neurotrophic tyrosine receptor kinase gene fusion and have been identified as an effective target for anticancer therapy. The study of the mechanism between polo-like kinase (PLKs) and pyrazol inhibitors was performed using 3D-QSAR modeling, molecular docking, and MD simulations in order to design high-activity inhibitors. The HQSAR (Q 2 = 0.793, R 2 = 0.917, R 2 ext = 0.961), CoMFA (Q 2 = 0.582, R 2 = 0.722, R 2 ext = 0.951), CoMSIA/SE (Q 2 = 0.603, R 2 = 0.801, R 2 ext = 0.849), and Topomer CoMFA (Q 2 = 0.726, R 2 = 0.992, R 2 ext = 0.717) showed good reliability and predictability. All models have been successfully tested by external validation, so all five established models are reliable. The analysis of the different contour maps of different models gives structural information to improve the inhibitory function. Molecular docking results show that the amino acids Met 592, GLU 590, LEU 657, VAL 524, and PHE 589 are the active sites of the tropomyosin receptor TRKs. The results obtained by MD showed that compound 19i could form a more stable complex protein (PDB id: 5KVT). Based on these results, we developed new compounds and their expected inhibitory activities. The results of physicochemical and ADME-Tox properties showed that the four proposed molecules are orally bioavailable, and they are not toxic in the Ames test. Thus, these results would provide modeling information that could help experimental researchers find TRK type I inhibitors more efficiently. Communicated by Ramaswamy H. Sarma


Introduction
In recent years, the most infectious and dangerous disease is cancer (Burnet et al., 1957), which is manifested by abnormal regulation of various pathways and remains the second most common disease and one of the major health problems in the world (Cancer, n.d.).Despite the development of materials and methods used to discover a drug, the problem of mortality from the disease is still increasing (Cancer Statistics, 2022, 2022).Among the methods used to discover drug candidates are computational methods that reduce the cost of drugs and increase their effectiveness.Receptor kinases of tropomyosin (TRKs), a family of three receptor tyrosine kinases (RTKs), regulate cell differentiation, proliferation, survival, and pain (Bernard-Gauthier et al., 2015).TRKs contain an extracellular domain for ligand binding, a transmembrane domain, and a short intercellular domain that consists of approximately 70 amino acids before and 15 amino acids after the kinase domain (Huang & Reichardt, 2001).TRKs have great potential in cancer drug research and development, and the design of TRKs with high efficiency and selectivity has become the research goal of drug workers (Drilon, 2019;Kciuk et al., 2022aKciuk et al., , 2022bKciuk et al., , 2022c)).At present, irreversible TRKs inhibitors have been introduced into clinical trials (Zhong et al., 2021).
A series of novel tetrahydropyrrolo [3, 4-c] pyrazol derivatives were synthesized and evaluated for their biological activity by Tianxiao Wu et al (Wu et al., 2021).This series is under molecular modeling study to propose further candidate TRKs inhibitors.Three-dimensional quantitative structure-activity relationship (3D-QSAR), molecular docking, molecular dynamics (MD), and (Chemical absorption, distribution, metabolism, excretion, and toxicity (ADMET) molecular property studies have been very important methods to generate predictive and robust models and also to identify the most important sites in the protein and verify the stability of each ligand, in order to predict new candidates for drug development.Hologram quantitative structure activity relationship (HQSAR), comparative molecular field analysis (CoMFA), comparative similarity indices analysis (CoMSIA), and Topomer CoMFA methodologies are used to correlate the dependent variations (pIC 50 ) with different molecular properties, such as steric and electrostatic chambers … ., in order to establish a good mathematical model (Keretsu et al., 2021).To verify the predictive ability and robustness of the different proposed models such as HQSAR, CoMFA, CoMSIA, and topomer CoMFA, we tested the obtained results by external and internal validation (Kovalishyn et al., 2021).We used molecular docking to predict the optimal binding configuration of a ligand and to understand the differences in structural interactions between a more active ligand and its protein target, and then we used MD simulations to analyze and further investigate the details of the interaction and stability of strong ligands in protein targets (Huang et al., 2019).
In order to explain the structure-activity relationship (SAR) and explore the possible optimization direction for such inhibitors, the 3-QSAR model was developed based on the docking poses.The contour maps could provide theoretical guidance for the design of new inhibitors TRKs.Based on the results of molecular docking and contour maps, four novel inhibitor molecules with high predictive activity were designed.Finally, we evaluated the ability of the four compounds to act successfully as drug candidates, tested by pharmacokinetic and pharmacodynamic parameters (ADMET and Lipinski rule).Thus, the present study will provide valuable insights for the design of new TRKs inhibitors.

Molecular modeling and alignment
The most crucial and sensitive stage for CoMFA and CoMSIA analysis is molecular alignment (Bringmann & Rummey, 2003).In this study, we created a straightforward alignment using the SYBYL-X 2.1 program.Each molecular structure created using the sketching module was reduced using the relevant Gasteiger-Huckel atomic partial charges on the SYBYL-X 2.1 platform while subject to the standard Tripos force field (Tsai et al., 2010).Additionally, the Powell gradient algorithm's convergence criteria were set at 0.005 kcal/mol, and 1000 iterations were required to reach a stable conformation (Sein et al., 1998).

Generation of 3D-QSAR models
Using Sybyl X-2.1 software (Tripos, Inc., USA), studies of the CoMFA and CoMSIA fields were carried out (Kang et al., 2004).The CoMFA descriptors were generated for each lattice with a grid spacing of 1 and extending to 4 units in all three dimensions within the given region (B€ ohm et al., 1999).Using the standard Tripos force field, Van Der Waals potential, coulombic terms, and an sp 3 hybridized carbon atom with a charge of þ1e as the probe atom, the steric and electrostatic fields' energies were calculated (Er-rajy et al., 2022).The fields' values were truncated at 30 kcal/mol (Sato et al., 2021).A distance-dependent Gaussian type of physicochemical characteristic was chosen for CoMSIA in order to prevent any singularities at the atomic position.Similar standard parameters, including those for steric (S), electrostatic (E), and hydrophobic (H) effects, were used for the CoMSIA field calculation, but no arbitrary constraints were created for the CoMFA field calculation (Sippl, 2002).Column filtering and attenuation factors, as well as attenuation factors, were both set to the default values of 0.3 and 2.0 kcal/mol, respectively (Guo et al., 2005).

Partial least squares analysis and validation
The electrostatic, steric, hydrophobic, hydrogen bond donor (D), and hydrogen bond acceptor (A) characteristics, each with polo-like kinase 1 (PLK1) inhibitory activity, were evaluated using the partial least squares (PLS) approach.The leave-one-out (LOO) cross-validation technique is employed in the PLS analysis to identify the ideal number of components (NOC), from which the final CoMFA and CoMSIA models were derived.We used the cross-validation correlation coefficient (Q 2 ), the coefficient of determination (R 2 ), the F-value (Fischer test), and the standard error of estimate (SEE) to determine the models' relevance (fadili, Er-Rajy, et al., 2022).Finally, using a non-validated PLS analysis, the various CoMFA and CoMSIA models were created.We require the biggest statistical Q 2 and R 2 values and the smallest SEE values to select a good module (Roy & Pratim Roy, 2009).
The biological activities of the external test set (8 chemicals) were predicted by the 3D-QSAR models in order to assess the predictive power of these models.

Dynamic molecular
Based on the molecular docking results, the two most active ligands were selected for MD simulations to determine the stability of the molecular bonds with the target protein (Kciuk et al., 2022dMujwar et al., 2021;Parihar et al., 2022;Shinu et al., 2022).The MD was run for 100 ns using the Desmond program through Maestro of the Schr€ odinger Suite, the OPLS3e force field (Miao et al., 2019).The OPLS3e force field was used for the MD simulation of two complexes and ligands bound (the most active) to the base receptor by solving it in a 10-size aqueous container (Kurki et al., 2022;Mujwar, 2021a;Mujwar & Kumar, 2020;Mujwar et al., 2022;Sharma et al., 2020).For a 10 region, the single point charge (SPC) model was used to infer crystallographic water molecules (TIP3P) under orthorhombic periodic boundary conditions (Bhatia et al., 2022;Bizzarri & Cannistraro, 2002;Singh et al., 2022).The use of the OPLS3e force field and the single point charge model for protein complexes using two small ligands (the most active) to obtain the best and most reproducible results (Chohan et al., 2016).
The system was neutralized by the addition of counter ions, including Na þ and Cl À , followed by the minimization of their energy (Pang, 1999).The long-range electrostatic connections between the protein and the complexed ligands were calculated using the Ewald particle mesh (PME) approach with a grid distance of 0.8 and a radius of 9.0 for coulombic interactions after running the MD simulation of the NPT package of Nose-Hoover thermostat and barostat.This was applied to keep the temperature (300 K) and pressure (1 bar) constant of the systems, respectively (Gupta et al., 2022a(Gupta et al., , 2022b;;Schulz et al., 2009).The specific binding interactions of the ligand with the viral protein were found using the interaction diagram module of the 2019-4 Desmond package simulation (Scarabelli et al., 2022).

Design new molecular, Lipinski's rule and ADMET prediction
From the results obtained by QSAR and docking methods, we proposed four new molecules.To calculate the different physicochemical properties and pharmacokinetic parameters of the new proposed compounds (Walters, 2012), we used the two websites, Swiss ADME prediction and pkCSM (PkCSM, n.d.; SwissADME, n.d.), Lipinski's five rules for oral bioavailability and ADME-Tox properties (Zhang & Wilkinson, 2007).

Alignments rules
To build a very strong and reliable model, you have to overlay all the molecules in a very sensitive way.All the molecule derivatives were added to a database.The database was superimposed on the common core by the rigid distillation method of Sybyl X-2.1, using compound number 19i (the very active one) as a model (Figure 1).

Results of HQSAR models
After the database superposition, we built ten HQSAR models by different combinations of the characteristic parameters, namely the number of atom fragments and the length of the hologram.We select the different models by the smallest error.Table 2 shows that the best HQSAR model is model number 5, obtained by combining two descriptors of atomic fragments (A), and bonds (B) by setting the default fragment length: minimum 4, and maximum 7, giving Q 2 cv and R 2 ncv , respectively, 0.739 and 0.945, with HL of 97 and N ¼ 6.
To improve the characteristics of the obtained HQSAR model, we change the fragment lengths in a regular way and observe the change in the characteristics of each model.
From the table, models 5-7 have the most significant statistical metrics by Q 2 cv and R 2 ncv , respectively, 0.793 and 0.917, with HL of 59, N ¼ 5, and fragment size from 6 to 9. The results of the biological activity predicted by the best obtained modulus (model 5-7) are presented in table 2.

HQSAR contribution maps analysis
To understand better the effect of each atom of the two selected molecules on the biological activity, the different colors of the diagram of the HQSAR model reflect an effect on the biological activity of each molecule.The green or yellow shows a positive contribution to the biological activity, and the orange and red colors indicate a negative contribution to the biological activity.The contribution diagram of compounds 19i (high activity) and 17t (low activity) is presented in Figure 2.
According to Figure 2, in molecule 17t, the radical morpholine and fluorine groups showed a negative contribution (red color) to the biological activity of the less active molecule due to electronegative groups such as the CF3 group substituent and the benzene group of morpholine substituents.
On the other hand, in molecule 19i, we observe that the morpholine group represents a positive effect on the activity (yellow color), and also in the main chain, we observe a positive effect on the nitrogen molecule (green color), and no effect on the biological activity.This gives an idea of the difference in biological activity between the two molecules.
To conclude, according to the analysis of the obtained results, we can say that to increase the inhibitory activity of each ligand, it is necessary to add donor groups and electropositive molecules.

CoMFA and CoMSIA studies
The table (Supplementary Table S2) gathers the observed and predicted activity and their residuals by the different models obtained (HQSAR, CoMFA, CoMSIA, and the Topomer CoMFA).The QSAR models were proposed to quantitatively explain and predict the anticancer activities of a series of tetrahydropyrrolo [3, 4-c] pyrazol derivatives.To determine the proper statistical parameters for each module (Olah et al., 2004), PLS cross-validation analysis was applied to the independent variables collected from the training set.The following table (Table 4) shows the results for the different statistical modules and their own statistical parameters.
In the CoMSIA study, the evaluation analysis of the seven selected models with different combinations of fields, such as stereoscopic (S), hydrophobic (H), electrostatic (E), and hydrogen bond donor (D) was performed.Table 3 shows that CoMSIA/SE and CoMSIA/SEHD were the best models among the various field combinations chosen, which obtained the highest Q 2 cv values of 0.603, and 0.634, respectively, with principal components 2, and 1, respectively, R 2 ncv values equal to 0.849, and 0.933, respectively, and a lower value of SEE equal to 0.312, and 0.344, respectively.To ensure the robustness of the proposed model, the results obtained are the objective of external validation.

External validation
The external validation allows us to evaluate and judge the predictive capacities of the obtained model; the following table groups the different external validation statistical parameters.
According to the results in Table 5, the five models HQSAR, CoMFA, CoMSIA/SE, CoMSIA/SEDH and Topomer CoMFA have better external prediction coefficients than the five models (R 2 pred Þ, which were used to validate the external predictive ability of the five models.The R 2 pred value of the five models is 0.793, 0.582, 0.603, 0.634 and 0.726, respectively.Therefore, the five models have Q 2 values greater than 0.5 (Golbraikh et al., 2003;Roy et al., 2012).
The external validation's predictive capability supported the five QSAR models' favorable stability estimates and high levels of predictive quality.We evaluated the best-chosen models' robustness and predictability in order to validate these findings (Er-rajy et al., 2022).Table 5 displays the outcomes of the five models' calculations for the external validation test.The results show that all five models met the Golbraikh and Tropsha requirements (Alexander et al., 2015),

CoMFA, Topomer CoMFA and CoMSIA contour map study
To see the details of the three 3D QSAR models contained as a starting point, we used the third and most active chemical in the sequence.The CoMSIA and CoMFA contour charts' steric, electrostatic, and hydrogen bond acceptor fields are shown in Figures 3 and 4, respectively.

CoMFA contour maps
The steric and electrostatic contour maps of the CoMFA model are shown in Figure 3.The CoMFA steric contour map (Figure 3 (Ster)) showed a combination of 80% green and 20% yellow contribution in the contour map with the most active compound 19i.
In Figure 3 (Ster), the green parts indicate that a steric contribution increases the potency of the molecule with a contribution of 80%, while the yellow parts indicate the steric hindrance that decreases the potency of the compound with a contribution of 20%.From the figure, we observe a small part of the green contours of the steric field surrounding the morpholino group, which is allows for an increase in the potency of the more active molecule.On the other hand, there is a large part and two small parts of yellow contours that surround the morpholino group, so steric hindrance that decreases the potency of the more active molecule surrounds the whole morpholino group.The electrostatic contour map of CoMFA is shown in Figure 3 (Electro).The red and blue contours represent 20% and 80%, respectively.A blue-colored region surrounded the 2,4-difluorobenzene group, ie the positively charged favorable 2,4-difluorobenzene group.The other part, a small red region located next to the methyl group is favorable to negatively charged groups.From the analysis of both steric and electrostatic contour maps, it can be concluded that the two morpholino groups and 2,4-difluorobenzene are responsible for the most active molecule.

CoMSIA/SEDH contour maps
The CoMSIA/SEDH contour maps are shown in Figure 4.The CoMSIA steric contour map (Figure 4 (Ster)) has the same  result as the CoMFA contour map, so the discussion of this contour map will be the same as the CoMFA contour map.
For the electrostatic CoMSIA contour map (Figure 4 (Electro)), there are two blue region clusters surrounding the 2,4-difluorobenzene group and the morpholino group, that is, the two positively charged groups, but no red region clusters.This means that the most active molecule does not have any negatively charged groups.We observed a small cyan contour in the middle of the molecule and another small orange contour in the contour map of the hydrogen bond donor (Figure 4 (H-b Donor)), indicating that the hydrogen bond donor has no significant influence on the inhibitory activity of the most active molecule.On the hydrophobic bond contour map (Figure 4 (Hydro)), the hydrophobic contour map of CoMSIA showed a combination of 80% purple (favorable) and 20% magenta (unfavorable) contributions in the contour map (Ai et al., 2010).From the figure, we observed two purple contours covering the two groupings of 2,4-difluorobenzene, and morpholino, which have a favorable effect on the inhibitory activity of the most active molecule.From the analysis of two contour maps, CoMFA and CoMSIA, we can say that the two groups morpholino and 2,4-difluorobenzene lead to an increase in the activity of the selected compound.

Topomer CoMFA
For the Topomer CoMFA, the model was generated by dividing the molecules into different fragments, aligning each fragment in a specific a manner using a model compound (compound 19i), and calculating the steric and electrostatic descriptors of the aligned fragments.
The 3D contour charts of the CoMFA topomer around position Ortho and Meta of cyclohexane were created by plotting the coefficients of the model.These maps are therefore presented using compound 19i as the reference structure.In the steric field of CoMFA topomer around position Ortho in cyclohexane, the green contours located at the nitrogen site (Figure 5(a)), and in the electrostatic field, the blue contours located at the nitrogen site indicate that the electropositive groups around the molecule would be favorable, and the red contours indicate that the electronegative groups would be favorable (Figure 5(b)).
In terms of fluorine group contours, yellow contours (Figure 5(c)) were located near fluorine, and green contours were absent, for red contours were located next to fluorine molecule and blue contours were located next to benzene ring meta position (Figure 5(d)); this demonstrated that a moderately bulky group with electronegative potential on fluorine site decrease the anticancer activity.

Molecular docking
To explain how the most active ligand will integrate with its protein, we performed a molecular docking of the two most active molecules with the protein (PDB id: 5KVT).In order to compare with the interactions established with the co-crystallized ligand (Entrectinib) shown in Figure 6, which indicates that amino acids Met 592, Glu 590, Leu 657, Val 524 and Phe 589 are the active sites of the target protein, as indicated by the ProteinsPlus online server.The results found by molecular docking for the two selected compounds (molecule 19i and 19k) are presented in Figure 7.
The first visual monitoring of two results shows that there is one hydrogen bond with residues ARG 654, with a distance equal to 2.75, and two halogen bonds with residues MET 592 and GLU 590, with a distance equal to 3.00 and 3.54, respectively, and two Pi-Anion bonds with residues ASP 668 and ASP 596, with a distance equal to 4.03 and 4.77, respectively, and one Pi-Sigma bond with residues PHE 589, with a distance equal to 3.47, and one Pi-alkyl bond, for the molecule 19i (Figure7(a)).In molecule 19k, we also observe three hydrogen bonds with the residues MET 592, ASR 596, and ARG 654 with a distance equal to 2.58, 2.84, and 2.35, respectively, a Pi-Sigma bond with the residue LEU 657, and a Pi-Alkyl bond (Figure 7(b)).Thus, the two molecules are attacked almost at the same site, which means that both molecules are inhibitors of the docked protein.We conclude that the amino acids Met 592, GLU 590, LEU 657, VAL 524, and PHE 589 are the active sites of the tropomyosin receptor protein kinases TRKs.

Dynamic simulation
The MD simulation is a useful technique to give an idea of the stability of the ligand in the active site of the protein  and the contribution of key amino acids in the ligand-protein interaction (Wang et al., 2018).Two of the most active molecules (19i and 19k) with the target protein were selected for MD up to 100 ns.The MD trajectories were examined using the RMSD parameter and the protein interactions with the ligand.
MD simulation has proven to be a powerful tool to study the stability of ligand docking postures and the contributions of key amino acids in proteins.In order to evaluate the stability of the complexes and to observe the possible ligand binding modes, we performed MD simulations of compounds 19i and 33 for 60 ns.As shown in Figure 8(a), the root mean square deviation (RMSD) value of the protein backbone showed a change in RMSD from 0.66 to 3.2 Å during the first 70 ns.Furthermore, the RMSD value of the protein configuration continued to increase during the simulation until it stabilized around 3.2 Å after 75 ns.In the case of the initially designed molecule in Figure 8(a), the RMSD ranged from 0.6 to 1.6 Å over 100 ns.The average RMSD of the crystal structure complex was about 3.2 Å, while that of the designed molecule was about 1.2 Å, clearly indicating that the designed molecule retains greater stability in the protein-ligand complex.For the second complex, as shown in Figure 8(b), the RMSD value of the protein backbone showed a change in RMSD from 1.00 to 2.5 Å throughout the simulation process, despite some small fluctuations authors of 78 ns.In the case of the designed molecule 19k Figure 8(b), the RMSD was between 0.6 and 2.7 Å during the 100 ns, around 20 ns there is a fluctuation going up to 2.7 Å, after there is another fluctuation around 40 ns going up to 2.4 Å and after becomes remains authors of 1.5 Å.The average RMSD of the complex crystalline structure was of about 2.1 Å, while that of the designed molecule was of about 1.3 Å.
Analysis of the contacts observed between the macromolecular complexes bound to the two most active compounds during a 100 ns MD simulation, according to Figure 9(a), revealed that the macromolecular structure had the amino acid Arg654 exhibited 90% hydrogen bonding interaction and 10% water bridge bonding, and the amino acid Asp674 exhibited 80% hydrogen bonding interaction and 20% water bridge bonding.For the second complex, the macromolecular structure had the amino acid Arg 654 exhibited 40% hydrogen bond interaction, 50% water bridge bond, and 10% hydrogen bond interaction.So, after the analysis of two figures, we conclude that the molecule 19i keeps greater stability in the protein-ligand complex.

Design of new compounds
The SAR information revealed by the above 3D-QSAR and HQSAR contour map analysis was summarized as shown in Figure 8, which may be useful for designing novel high activity TRKs inhibitors.
From the results obtained by a different method (Figure 10) and the discussion we gave, based on the SAR information obtained by the different contour maps CoMFA and CoMSIA to predict new designed compounds (Table 6).The new compounds were included in the test set, and their activities were predicted by CoMFA, CoMSIA, and HQSAR analyses.
After the changes in the groups that influence the inhibitory activity of the proposed molecules, we can calculate the predicted activity of four proposed compounds by the obtained equations (CoMFA, CoMSIA, and HQSAR).The obtained results of the predicted activities of four new compounds (1, 2, 3, and 4) are presented in Table 6.From Table 6, all four molecules have predicted activity (pIC 50 ) greater than 7 for the three proposed QSAR models, thus the four proposed molecules have greater inhibitory activity.
Then for more information about the different pharmacophore present in the most active molecule in the proposed compound, we selected the generated models such as CoMFA and CoMSIA to visualize the different pharmacophores of the selected compound.
The pharmacophore of the most active molecules (molecule 4), in proposed molecules are presented in two following figures.
In Figure 11 (left), we observe a large part of green contours of the steric field surrounding the morpholino group, which allows an increase in the potency of new most active molecule.On the other hand, there is a large part and two small parts of yellow contours surrounding the -BH2 group, so the steric hindrance that decreases the potency of new more active molecule surrounds the whole morpholino group.The electrostatic contour map of CoMFA is shown in Figure 11 (right).A blue colored region surrounds the -BH2 group, which we added that increases the strength of new molecule.The other part, a small red region located next to the -BH2 and morpholine is favorable to negatively charged groups.The analysis of the steric and electrostatic contour maps allows us to conclude that the two morpholino groups and -BH2 are responsible for the strength of new molecule.
From the CoMSIA-SEDH contour (Figure 12), we can say that the different moieties such as morpholino and -BH2 lead to an increase in the activity of the newly selected compound.
Therefore, after the analysis of different map contours of the most active molecule and also for the new molecule, it is observed that the two substituents morpholine, fluorine and -BH2 are responsible for the strength of the activities of two molecules.

Lipinski's rule and in silico ADME-Tox evaluations
To verify Lipinski's rule of similarity properties of drugs and ADME-Tox properties, we used two websites, such as Swissadme and pkCSM, to calculate the different properties of each molecule we proposed.
According to the results obtained in Figure 13 and Table 7; Figure 13 shows that all the newly proposed compounds have oral bioavailability properties in the pink zone, which represents the optimal range for each property, and Table 7 affirms this result.In conclusion, we can say that four compounds have a low attrition rate during clinical trials.To avoid falling into clinical trials because of their toxicity or poor pharmacokinetics, after checking for similarity to the drug, they were subjected to ADMET prediction.The results obtained are presented in the following table.From Table 8, it can be observed that all the water solubility values are less than zero, so it can be said that the four new compounds are highly soluble in water.The gastrointestinal absorption capacities of the four compounds were strong; compounds 1, 2, and 4 were unable to cross the BBB, but compound 3 was able to cross the BBB.All the compounds were glycoprotein substrates, and all four compounds could not inhibit the cytochrome P450 subtype 1A2 enzyme, but all the compounds inhibited the cytochrome P450 subtype 2D6 enzyme.All four compounds tested were cytochrome P450 enzyme substrates.The four novel chemicals are not dangerous, contrary to routinely used tests for toxicity such as the Ames toxicity test and the skin sensitization test.
In conclusion, based on the drug similarity studies and ADMET, we propose compound 4 as a TRKs enzyme inhibitor because it almost checks all the similarity properties of a drug, but it is necessary to conduct further studies for the area where this compound is a drug.

Conclusion
The study of the mechanism between PLKs and pyrazol inhibitors was performed using 3D-QSAR modeling, molecular docking, and MD simulations in order to design high-activity inhibitors.The HQSAR, CoMFA, CoMSIA/SE, CoMSIA/SEHD, and Topomer CoMFA modules showed good reliability and predictability.The five established models are evaluated by external and internal validation methods with success.The analysis of the different contour maps of different models gives structural information to improve the inhibitory function.In addition, the interactions between the two selected compounds and the target protein were modeled by molecular docking studies, and the MD simulations up to 100 ns show the stability and provide information about the ligand-receptor interactions.Based on these results, we developed four new compounds were designed in silico, and we predicted their expected inhibitory activities by different model, and also their pharmacophore properties.The results of physicochemical properties and ADME-Tox showed that the four proposed molecules are orally bioavailable; furthermore, they are not toxic in Ames test.The results indicate that compound 4, present a good potential to become TRKs inhibitors with a better activity than the most active compound in the different molecules studied.Therefore, these results would provide the necessary physicochemical and pharmacokinetic properties and relevant information in drug discovery analysis.

Disclosure statement
No potential conflict of interest was reported by the authors.The ADME prediction results for the four new compounds, including water solubility, gastrointestinal absorption, blood-brain barrier (BBB) permeability, P-glycoprotein substrate, inhibition capacity, and cytochrome P450 enzyme substrate, obtained from the Swiss ADME prediction are listed in Table 8.Then, the pkCSM website was used to calculate the ADME/Tox properties.

Figure 1 .
Figure 1.Super-positioning of the core structure and alignment of the training and test set.
c Cross-validated correlation coefficient.d Cross-validated and non-cross-validated standard error.e Non-cross-validated correlation coefficient.

Figure 5 .
Figure 5. CoMFA Topomer model contour maps of the most active compound.

Figure 7 .
Figure 7. Two and three-dimensional binding interactions of the most potent compound 19i (a) and 19k (b) against TRKA kinase active site.

Figure 6 .
Figure 6.Two and three-dimensional binding interactions of Entrectinib against TRKA kinase active site.

Figure 8 .
Figure 8. RMSD observed during MD simulation of 100 ns for the macromolecular complex of compound 19i (a) and 19k (b) against human TRKA kinase enzyme.

Figure 9 .
Figure 9. Binding contacts observed during the MD simulation for 100 ns between the two best compounds 19i (a) and 19k (b) against human TRKA kinase enzyme.

Figure 13 .
Figure 13.Bioavailability radar charts of now compounds proposed.
Abbreviations: ID: Number of molecular; MW: Molecular weights; N Rot-bonds: Number of rotation bonds; N-HBA and N-HBA: Number of hydrogens bond acceptors and donors; TPSA: Topological polar surface area; Log P: Partition coefficient; S.A: Synthetic accessibility; Log S: Stands for solubility.

Table 1 .
Derivatives of benzylidenecyclohexenone and biological properties.Homepage, n.d.).The active site is surrounded by the three-dimensional configuration (X ¼ 27.578, Y ¼ À 29.337, Z ¼ 1.503), which has a grid with a size of 60 � 60 � 60 points in the x, y, and z directions.Following the production of the ligand and receptor, we performed molecular docking on two of the most active compounds

Table 2 .
HQSAR statistical metrics for different fragment distinction.

Table 3 .
HQSAR statistical metrics for a diversity of fragment sizes.
which means that they have improved predictive ability for new chemicals and comply with all validation techniques.

Table 5 .
Recapitulation of some statistical parameters.

Table 6 .
Design of new molecules and activity predict.

Table 7 .
The result of Lipinski's rule calculation for now compounds.

Table 8 .
The ADME prediction results for now compounds.ID Water solubility Gastrointestinal absorption BBB permeant P-glycoprotein substrate