^tButyl 4-(6-chloropyrimidin-4-yl)piperazine-1-carboxylate 1 was prepared via the reaction of 4,6-dichloro pyrimidine and Boc-protected piperazine following the reported method (29). In a 50 mL sealed tube containing mixture of dioxane/H₂O (8/2), a stoichiometric amount of 1 (1.0 mmol) and different arylboronic acid (2, 1.5 equiv.) was added followed by the addition of base (K₂CO₃, 6 equiv.). The reaction mixture was purged with nitrogen gas for 2 min, prior adding the Pd(0) catalyst (5 mol%) under a nitrogen atmosphere. The reaction was allowed to be heated at 100°C for 3 h. Upon the completion of the reaction, the solvent was evaporated on a rotary evaporator and the crude product was purified using gravity column chromatography (eluent: EtOAc/Hex = 8/2) yielding ^tbutyl 4-(6-arylpyrimidin-4-yl)piperazine-1-carboxylate 3. TFA-induced deprotection in dichloromethane (DCM) yielded compound 4-aryl-6-(piperazin-1-yl)pyrimidine 4.

To 50 mL of the round bottom flask containing 10 mL of DCM, compound 4 (1.0 mmol) was added while stirring. This is followed by portion-wise addition of 1H-pyrazole-1-carboximidamide (1.5 equiv.) and Et₃N (1.2 equiv.) at 0°C. The reaction was stirred at room temperature for 3−12 h. Following the completion of the reaction, the reaction mixture was diluted with diethyl ether. The resulting precipitate was filtered and washed with several portions of diethyl ether and dried to obtain 4-(6-(arylpyrimidin-4-yl)piperazine-1-carboximidamides 5−10 in good to excellent yield.

Purified MARK4 was used to conduct ATPase inhibition experiment following the reported method (18, 30). The concentration of the purified protein used was 5 μM, the ligand concentration was 0–10 μM and the substrate (ATP) was 200 μM. The ligand concentration was increased gradually against a given concentration of MARK4. The experiment involved combining MARK4 with freshly prepared ATP, resulting in a final reaction mixture volume of 100 μL. This mixture was incubated for 1 h at 25°C. This is followed by the addition of Malachite green solution (200 μL) to quench the reaction. The samples were allowed to incubate at room temperature for color development. Finally, the reaction mixture (100 μL) was filled into 96-well plates and analyzed using a multiscan ELISA reader (λ = 620 nm). All experimental reactions were conducted in triplicate to ensure statistical reliability.

A fluorescence-based binding assay was conducted to evaluate the binding efficiency of ligands. During the titration, the concentration of the pure MARK4 protein was kept constant while increasing concentrations of the ligand (31). The emission data were collected using a Jasco FP-6200 spectrofluorimeter (Japan), with a slit width of 10 nm and medium sensitivity. The samples were excited at 280 nm, and the emission spectra were recorded. Blank subtracted spectra are reported, and each measurement was performed in triplicate. For fluorescence values, we took the inner filter effect into account (32). The binding constant (K) was calculated using the Stern-Volmer (SV) (Equations 1, 2).

Where F₀ denotes the fluorescence intensity of the protein in the absence of ligands, whereas F indicates the fluorescence intensity when ligands are present. K represents the binding constant, n indicates the number of binding sites, and C refers to the concentration of the ligand.

The MARK4 protein (PDB ID: 5ES1) was chosen as the target in this study (33). The 3D structure of MARK4 in complex with the ligand N-(6-amino-2,2-difluorocyclohexyl)-5-ethyl-4-(6-(trifluoromethyl)pyrazolo[1,5-a]pyrimidin-3-yl)thiophene-2-carboxamide (5RC) was retrieved from the Protein Data Bank (PDB) and saved in .pdb format.

MD simulation is a powerful technique used to study protein-ligand complexes, providing detailed insights into their dynamic behavior, structural stability, and interactions under physiological conditions. We utilized the Desmond module of Schrödinger for the simulations (41, 42). Using the system builder module, the ligand-protein complex was first selected and prepared for the MD simulation. An orthorhombic simulation box of dimensions 10 × 10 × 10 Å was created, and the SPC (Simple Point Charge) model was chosen as the solvent model (43). The OPLS3 force field was employed to represent the interactions within the system (44, 45). To mimic physiological conditions, sodium (Na⁺) and chloride (Cl⁻) ions were added to the system (0.15 M). Following this, energy minimization was performed for 100 picoseconds (ps) to relax the model into a local energy minimum. The minimized structure was then subjected to MD simulation for 500 nanoseconds (ns) under an NPT ensemble, which maintains constant temperature (300 K) and pressure (1.01325 bar) maintained using the Nosé-Hoover thermostat and Martyna-Tobias-Klein barostat, as implemented in the Desmond module (46–48). The results of the simulation were analyzed by comparing root mean square deviation (RMSD), root mean square fluctuations (RMSF) and secondary structure elements (SSEs).

The drug-likeness, absorption and distribution and toxicity of the compounds under investigation were predicted by SwissADME (49), pKCSM (50) and Qikprop (51, 52) tools.

Geometry optimization, chemical reactivity descriptors, frontiers molecular orbitals, energy and other properties of the compounds were computed by DFT (B3LYP/6-31G*) using Spartan 24 software (53).

Six-membered pyrimidines are a privileged scaffold in drug discovery, especially in neurodegenerative diseases (56). While 4,6-disubstituted pyrimidines have been extensively studied in the past against AD-related targets such as AchE (57), tau (58), β-amyloid aggregation (59) and monoamine oxidase (60), limited studies have been carried out against MARK enzymes. We recently demonstrated that the MARK4 inhibitory activity depended on the substitution pattern of 4-(4-(arylsulfonyl)piperazin-1-yl) flank in 4,6-disubstituted pyrimidine derivatives (29). To assess the inhibitory potentials of compounds 5–10, an ATPase inhibition assay against MARK4 was carried out following the reported method (18, 30). MARK4, currently does not have an FDA-approved inhibitor. However, numerous compounds are progressing through the pre-clinical stages of development. For this study, we chose Rosmarinic acid (IC₅₀ ~ 6.20 μM) as a reference inhibitor (25). Its potent bioactivity makes it a valuable reference point in exploration of potential MARK4 inhibitors. The results of the kinase activity assay are compiled in Table 1 and shown in Figures 4a–c. Regarding the activity, compound 5 (IC₅₀ = 5.35 ± 0.22 μM) showed the best activity, followed by 9 (IC₅₀ = 6.68 ± 0.80 μM). Notably, neither of these top-performing compounds contains any halogen substituents. In contrast, compounds 6, 7, and 8, which feature two or more halogen atoms attached to their aromatic rings, displayed lower activity, with IC₅₀ values of 11.76 ± 1.19 μM, 16.53 ± 1.71 μM, and 15.78 ± 0.92 μM, respectively. Furthermore, compound 10, which possesses a trifluoromethoxy group at the 4-position, showed intermediate activity with an IC₅₀ = 7.28 ± 0.49 μM, situating it between the two extremes observed. To further elucidate the underlying reasons for these activity trends, we performed extensive in-silico studies (vide-infra). Following this, we evaluated the effects of compounds 5 and 9 across a range of concentrations (0 to 10 μM). Notably, we observed a dose-dependent activity pattern: as the concentration of these compounds increased, there was a corresponding decrease in their activity. This reduction in MARK4 activity is noteworthy and suggests that the behavior of compounds 5 and 9 is comparable to, or even superior to, that of several small molecules derived from both natural sources (25, 26) or synthetic origin (14). On the contrary, when we conducted assays targeting pyruvate dehydrogenase kinase 3 (PDK3), these compounds displayed a weaker affinity (data not included). Such selectivity towards MARK4 of the azole-based compounds is already known (16). Therefore, compounds 5 and 9 stand out as representatives of a novel class of MARK4 inhibitors.

The presence of aromatic amino acid (tryptophan) in the MARK4 protein results in a characteristic fluorescence emission near 344 nm when excited at 280 nm (25). To get deeper insight into the interactions between the 5-MARK4 and 9-MARK4 complexes, we conducted fluorescence binding studies and determined the binding constant (K). As we gradually introduced increasing concentrations of compound 5 (0–7 μM) to the native MARK4, we observed a notable decrease in the fluorescence intensity (Figures 5a,b). This shift was accompanied by a slight blue shift in the emission maximum, approximately 5–7 nm. A similar pattern emerged with the 9-MARK4 complex, further supporting the notion that the ligand is effectively binding to the receptor. A decrease in intensity is also indicative of a stable complex formed between the ligand and protein (61). By employing a modified Stern-Volmer equation, we calculated the binding constants (Figures 5c,d), which were found to be 1.5 ± 0.51 × 10⁵ M⁻¹ for 5 and 1.14 ± 0.26 × 10⁵ M⁻¹. These binding constants highlight a strong affinity between the compounds and the MARK4 protein.

We conducted molecular docking studies to underpin the role of different functionalities on their biological activity and elucidate the binding mode and strength between MARK4 and ligands. The ligand (5RC) was used as the reference compound. Table 2 collects the docking score and binding residues of MARK4 complexed with 5RC and compounds 5–10. All compounds, including the inhibitor 5RC, are targeted in the ATP binding pocket and bind to the deep binding cavity of MARK4 (62). The docking score of the compounds was higher than the ligands like huperzine (19). The residues responsible for binding the ligands to the active site of MARK have been reported to be involved in the binding of ATP (a pyrimidine-based molecule) and other inhibitors (62). Figures 6, 7 displays the 2D/3D docking poses and interacting residues between 5RC or 5 and the MARK 4 receptor, while the other compounds are given in Supplementary Figures S9, S10. 5RC formed H-bond with Lys88, Ala138, Glu185, Asn186, and Asp199 of the active site of MARK4 (Figure 6). Binding free energy (ΔG) for MARK4 complexed with 5RC was −8.0 kcal/mol, indicating reasonable binding affinity. However, an increase in affinity (ΔG = −8.4 kcal/mol) can be noted for the complex containing 5. Indeed, the NH₂ group of guanidine facilitated the formation of H-bonds with the backbone amino acids Glu142 and Asp145 (Figure 7). Such residues have also been reported to be involved in other heterocyclic compounds (62). In addition, salt bridge interaction was also observed between the ligand and the side chains (Glu142 and Asp145), enhancing the ligand’s affinity for the kinase domain. A decrease in affinity (ΔG = −7.7 kcal/mol) can be seen when the tolyl ring attached to the pyrimidine was replaced by 3,4-dichlorophenyl in 6. It formed one only H-bond with Ser139 of MARK4. However, when the chloro was replaced with F (compound 7), an increase in the docking score (ΔG = −8.2 kcal/mol) was noted despite having one H-bond with Asp142 residue. However, the omission of a second chlorine and the introduction of a 4-trifluoromethyl (compound 8) or 4-trifluoromethoxy (compound 10) were found to be detrimental, as indicated by reduced docking score (ΔG = −7.8 kcal/mol for both), despite having one H-bond with Tyr137 and Lys67, respectively. Compound 9, containing an electron-rich benzothienyl core, demonstrated affinity like compound 8 (ΔG = −7.9 kcal/mol) and formed one H-bond with residue Tyr137. Summarily, this study found that among the screened compounds, 5 is the most promising candidate. It is also notable that the guanidine moiety attached to piperazine participated in H-bonding. A change in the aryl ring only affected the orientation of ligands of the molecules.

The MD simulations of native and MARK4 complex with compound 5 were conducted for 500 ns each. Both systems achieved equilibrium within the first 50 ns. The initial trajectories up to 150 ns displayed only minor fluctuations. However, after 150 ns, the native MARK4 exhibited a progressively increasing RMSD, indicating a deviated and less stable trajectory. In contrast, the 5-MARK4 complex demonstrated stable behavior beyond 150 ns, with the trajectory showing minimal fluctuations, suggesting that compound 5 stabilizes the MARK4 structure during the simulation. The average RMSD values for the entire simulation were 2.64 Å for the native, 2.35 Å for the 5-MARK4 complex, and 2.30 Å for the cocrystal ligand 5RC (Figure 8a). This further supports the conclusion that compound 5 enhances the structural stability of MARK4. To assess the flexibility of individual amino acid residues during the simulation, RMSF plots of the native, 5-MARK4 complex and cocrystal ligand 5RC were analyzed (Figure 8b). A significant fluctuation can be seen in the residues Gly66-Lys72, Ser210-Cys216, Gln228-Gly234 of the native MARK4, leading to RMSF values reaching ~6 Å with an average RMSF of 1.46 Å. The residue between Lys67 to Ala71 which forms the roof of the catalytic cleft, is already known for its flexibility (62). When complexed with compound 5, the fluctuation across these residues reduced, indicating the stability of MARK4 in the presence of the ligand. A low average RMSF (1.14 Å) is an indication that the protein was stable throughout the MD simulation. Similarly, the average RMSF of (1.13 Å) for the cocrystal ligand 5RC is almost like the compound 5.

Following the completion of the RMSD and RMSF analyses, we assessed the secondary structure elements (SSE), with particular emphasis on α-helices and β-strands (Figure 9). In the case of the native protein, 31.97% of the residues formed α-helix while 13.11% formed β-strands with a total of SSE content 45.08%. A slight drop in α-helix (31.53%) with an increase in β-strand (14.09%) and total SSE (45.61%) was noted for the 5-MARK4 complex. Similarly, the 5RC complex (31.70%) of the residues formed α-helix while (13.64%) formed β-strands with a total SSE content of 45.34%. This indicates that while the secondary structure of MARK4 remains largely conserved in the presence of compound 5, there is a slight increase in β-strand content, indicating a subtle stabilization of the secondary structure during the MD simulation.

A widely accepted methodology for predicting the potential success of new drug candidates involves the evaluation of their drug-likeness, bioavailability, and toxicity through the use of in-silico tools (63). This evaluation also helps determine how promising a drug candidate may be for further development and eventual clinical use. In literature, various tools have been reported that are used to estimate these properties (64, 65). In the present work, we computed these features using SwissADME (49), Qikprop (51, 52), and pKCSM (50) tools (Tables 3–5).

The results of the drug-likeness of the compounds 5–10, including the reference 5RC, computed using SwissADME (49) is given in Table 3. As can be seen, all compounds displayed excellent drug-likeness. None of them violated the parameters required for Ro5/Ro3 (66). Interestingly compound 5, which also displayed the best activity in vitro, showed the lowest topological polar surface area (TPSA = 82.1; limit ≤140 Ǻ²), molecular weight (MW = 296.4 amu; limit ≤500) and logP (1.38; limit logP ≤ 5). Notably, the predicted properties were better for the 5–10 than the reference compound 5RC. RADAR plot, in which each vertex shows optimum values of different descriptors, also did not show any deviation (Supplementary Figure S11). Among the main results of QikProp, all compounds show reasonable CNS activity (−1). The average number of HBD and HBA was found to be 3 and 4.5, respectively, which is within the recommended value (67). Solubility profiles like QPlogS (−4.79 to −5.69) and CIQP logS (−4.12 to −5.48) were excellent for all the compounds (recommended value: between −0.5 to −6.5). Similarly, the predicted brain/blood partition coefficient (QPlogBB) and apparent gut-blood barrier permeability QPPCaco (463.10 to 484.04) was high for the studied compounds. The most active compounds 5 and 9 showed low dipole moment and high human oral absorption (%HOA).

Despite having these advantages and drug-like features, we also noted some predicted properties that could undermine their potential. For example, the EGG plot (Supplementary Figure S12) indicates that all compounds are predicted to be well-absorbed by the gastrointestinal tract but have no access to the brain. Besides, PGP is also predicted to pump them out. We computed the toxicity profile of the compounds (Supplementary Table S1). All compounds are predicted to exhibit AMES toxicity, hERG-II toxicity and possible hepato-toxicity (50). QikProp also predicted IC₅₀ values for blockage of HERG K⁺ channels. These issues are concerning and will be studied in future investigations.

Density functional theory (DFT) calculations, like molecular docking and simulation studies, have become routinely used tools in drug discovery programs (68). A plethora of research available shows that DFT can be used in drug modelling (63), designing polymer based drug delivery (69), and calculating numerous other molecular and quantum descriptors (63). In the present study, we computed the energy and distribution of frontiers molecular orbitals, which was then used to estimate the molecular chemical indexes like chemical reactivity, hardness, softness, etc. for the synthesized 5–10 and reference (5RC) compounds. In 5RC, highest occupied molecular orbital (HOMO) was located over 6-amino-2,2-difluorocyclohexyl unit while the lowest unoccupied molecular orbital (LUMO) was over 6-(trifluoromethyl) pyrazolo[1,5-a]pyrimidin-3-yl fragment (Figure 10). The HOMO-LUMO gap, which is a key indicator of a molecule’s chemical stability and reactivity, was found to be 3.82 eV (Supplementary Table S2). In the case of compounds 5–10, both HOMOs/LUMOs were found to increase. Such band gaps have been reported in small molecules recently (70). We can see that the HOMO level of the reference compound (5RC) was deeper, while LUMO was higher than the compounds 5–10. It is also interesting to note that the most active compounds 5 (−5.71 eV) and 9 (−5.74 eV) have the highest HOMO levels. Among the studied compounds, (5) exhibited the highest HOMO-LUMO gap (4.44 eV). Like the reference compound 5RC, here, too, HOMO was located over piperazine-1-carboximidamide while LUMO over 4-(6-(p-tolyl)pyrimidin-4-yl) fragment. Similar observations were noted for the other compounds of series (Supplementary Figures S13, S14). Red regions over the molecules in the electrostatic potential (ESP) map indicate electron-rich region, which participated in the H-bonding, as revealed in the docking studies.