Computational screening of natural inhibitors against Plasmodium falciparum kinases: Toward novel antimalarial therapies

TL;DR

This study uses computational methods to find natural compounds that could act as new antimalarial drugs by targeting key proteins in Plasmodium falciparum.

Contribution

The study introduces a novel computational screening approach to identify natural inhibitors against P. falciparum kinases for potential antimalarial therapies.

Findings

Ligand-13 showed strong binding to CK2 with a score of −11.468 kcal/mol and demonstrated good pharmacological and toxicological properties.

Ligand-9 was identified as a strong inhibitor of PKG with a score of −7.490 kcal/mol.

Density functional theory and molecular dynamics simulations confirmed Ligand-13's stability and reactivity.

Abstract

An important worldwide problem is the resistance of Plasmodium falciparum to practically all antimalarial medications. Therefore, new treatment approaches are urgently needed. The development of antimalarial medications frequently involves two important therapeutic targets: casein kinase 2 (CK2) and cGMP-dependent protein kinase (PKG). To identify naturally occurring chemicals that could be used as antimalarial medications to combat multidrug-resistant P. falciparum, we used a multi-targeted in silico strategy in this study. The top 20 compounds, including the reference drug RY-1–65, were selected after pharmacophore-based virtual screening of naturally produced compounds. These compounds were subsequently docked onto both target proteins using Maestro (Schrödinger 2020−3). The best-scoring compounds against PKG and CK2 were Ligand-9 (−7.490 kcal/mol) and Ligand-13 (−11.468 kcal/mol), respectively. These lead compounds may be useful as therapeutic targets based on an assessment of their pharmacological, toxicological, and bioactivity characteristics. Furthermore, Ligand-13’s strong reactivity and stability were demonstrated by density functional theory analysis, and these findings were confirmed by molecular dynamics simulations and binding free energy MMGBSA calculations. These results imply that Ligand-13 may be a promising antimalarial medication.

Figures

The workflow of the in silico-based methodology.

Selected target proteins of Plasmodium falciparum.A: cGMP-dependent protein kinase PKG. B: Protein kinase CK2 catalytic domain.

The 3D structures of the reference compound (U) and top 20 pharmacophore-based virtually screened naturally derived compounds (Ligand-1 to Ligand-20 presented chronologically as A to T).

Visual representation of the site binding of Target-1 (PKG) to the reference drug and Ligand-9.A: 3D representation of PKG protein. B: PKG–Ligand-9 3D and 2D interaction with labeled interacting residues. C: PKG–reference drug 3D and 2D interaction with labeled interacting residues.

Visual representation of the site binding of Target-1 (CK2) to the reference drug and Ligand-13.A: 3D representation of CK2 protein. B: CK2–Ligand-13 3D and 2D interaction with labeled interacting residues. C: CK2–reference drug 3D and 2D interaction with labeled interacting residues.

Boiled egg diagram of the reference drug and Ligand-9 and Ligand-13, in which molecule-1 refers to the reference drug, molecule-2 refers to Ligand-9, and molecule-3 represents Ligand-13.The radar charts in the left corner show A: reference drug; B: Ligand-9; and C: Ligand-13.

Bioactivity of the reference drug (A), Ligand-9 (B), and Ligand-13 (C).

Visual representation of the molecular electrostatic potential and HOMO–LUMO energy gap of the reference drug (A), Ligand-9 (B), and Ligand-13 (C).

Root mean square deviation (RMSD) of the docked complex of CK2 with reference drug (A) and with Ligand-13 (C), with the highlighted superimposed structures B and D, respectively.The RMSD plot of the protein–ligand complex displayed two notable fluctuations around 16 ns and 54 ns, indicating transitions between distinct conformational states. Based on the RMSD trajectory, the system remained stable during 0–16 ns (Basin I), shifted to a new conformation from 16–54 ns (Basin II), and then returned toward the initial stable state after 54 ns (Basin I).

A, B, and C plots of the docked complex of CK2 with Ligand-13 (blue) and the reference drug (purple) show the root mean square fluctuation (RMSF), radius of gyration (Rg), and solvent-accessible surface area (SASA), respectively.

Comparative representation of principal component analysis (PCA), free energy landscape, and dynamic cross-correlation of the docked complex of CK2 with the reference drug (A, C, E) and with Ligand-13 (B, D, F).

Comparative visualization of the directional movements of protein during the docking of CK2 with the reference drug (A, B) and with Ligand-13 (C, D) observed at different intervals.The detailed reference structure of the protein is presented in the center.

Comparative analysis of ΔE and ΔG values of the reference drug (purple) and Ligand-13 (blue).

Equations

Full text

1. Introduction

An estimated 249 million clinical cases of malaria and 608,000 deaths due to the disease were reported globally in 2022. The World Health Organization’s most recent World Malaria 2024 report estimates that there were 597 000 malaria fatalities and 263 million cases globally in 2023. This translates to almost the same number of deaths and approximately 11 million more cases in 2023 than in 2022. In the WHO African Region, where many at-risk individuals still do not have access to the services necessary to prevent, detect, and treat the disease, approximately 95% of the deaths occurred [1]. Malaria is a vector-borne disease caused by a protozoan organism in the genus Plasmodium [2]. P. falciparum is by far the most virulent as well as the most prevalent species worldwide [3]. It is recognized as the deadliest malarial species due to the high number of deaths associated with it [4]. In most malaria-endemic nations, the World Health Organization recommends artemisinin-based combination therapy (ACT) as the first-line treatment for uncomplicated malaria [5,6]. The control of malaria is seriously threatened by the rise of antimalarial medication resistance in P. falciparum. Artemisinin resistance has developed and is mostly found in Southeast Asia, with few instances in Africa [7]. It is interesting to note that this drug resistance has spread to include other medications taken in combination with ACT. Moreover, resistance to commonly used treatment medications has also resulted from multidrug-resistant Plasmodium strains [8]. Malaria’s prevalence and drug resistance highlight the need to address enduring treatment hurdles despite efforts to eradicate the disease.

Multi-targeted therapy is necessary to eradicate P. falciparum because of its intricate life cycle [9]. Casein kinase 2 (CK2) and cGMP-dependent protein kinase (PKG) are two important kinase proteins used for this purpose. In P. falciparum, PKG is the primary effector of cGMP signaling. Additionally, in the erythrocytic and pre-erythrocytic embryonic stages of P. falciparum, it controls both the sexual and asexual cycles [10]. It has also been recognized as a possible target for antimalarial drugs in earlier studies [11,12]. CK2 in P. falciparum comprises three chains: two β-chains (B and C) for regulatory roles and one α-chain (A) for catalytic activity. The nucleus contains the alpha chain of CK2 (CK2α), which is specifically implicated in chromatin phosphorylation [13]. Furthermore, CK2 plays a role in the development of both the sexual and asexual blood stages [14].

The identification and preliminary characterization of a novel imidazole-based chemotype, RY-1–165, have been documented in previous studies. This chemotype exhibits promising cellular activity, good in vitro PfPKG inhibition, a correlation between in vitro enzymatic activity and efficacy, no structural alerts linked to genotoxicity, and no hERG problems, similar to other chemotypes [15]. The cyclopropyl group of compound 14a (RY-1–165) resulted in a three-fold increase in potency (IC_50_ = 100 nM, E_50_ = 7.6 μM) as a PfPKG inhibitor.[16].

The antimalarial drug of the future should ideally be a single-dose, cost-effective therapy capable of addressing parasite resistance. These conditions could be fulfilled through novel, easily synthesized natural compounds or their derivatives with selective antimalarial activity. Natural products such as chloroquine from the South American cinchona trees and artemisinin from the Chinese Artemisia annua are effective antimalarial drugs [17,18]. Despite resistance to these drugs, there is a need for further exploration of natural biotherapeutics.

In silico methods, which save time and money while enabling the assessment of novel molecule properties, such as toxicity and efficiency, before synthesis, can accelerate drug development [19]. Molecular docking and molecular dynamics (MD) simulations are two essential methods in molecular biology and computer-aided drug design. Molecular docking is frequently used to predict the binding orientation of therapeutic candidates to their protein targets to estimate the affinity and activity of small molecules [20]. In P. falciparum, PKG is the primary effector of cGMP signaling. Previous studies have identified PKG as a target, shedding light on the potential of inhibitors to prevent malaria. For instance, the interaction and stability of PKG-inhibitor complexes have been examined using docking and molecular dynamics simulations [21].

Therefore, this study employed an extensive in silico approach to explore and identify novel naturally derived compounds that can act as inhibitors against P. falciparum’s PKG and CK2 proteins to selectively target the malarial life cycle to safely and cost-effectively overcome antimalarial resistance.

2. Methodology

The workflow of this in silico-based study is outlined in Fig 1.

The workflow of the in silico-based methodology.

2.1. Target proteins

The target proteins were retrieved in Protein Data Bank (PDB) format from protein data bank (PDB) database [22]. Target-1, PKG (PDB ID: 8EM8), consisted of 857 amino acids and a molecular weight of 98.63 kDa [16]. Target-2, CK2 (PDB ID: 5XVU), consisted of 996 residues, a molecular weight of 121.58 kDa, and 3 chains (A, B, C) [23].

2.2. Reference inhibitor

The reference inhibitor for pharmacophore-based virtual screening (PBVS) of naturally derived compounds for both target proteins was RY-1–65 (PDB ID: 8EM8), a trisubstituted derivative of the antimalarial drug imidazole [16,15]. It was retrieved via the PubChem (CID: 164628676) database in SDF format.

2.3. Pharmacophore-based virtual screening

The PBVS was performed via Pharmit [24] using the crystal structures of the target proteins and the selected reference inhibitor. Pharmit ZINC Natural library, consisting of 1,441,177 conformers of 116,363 molecules, was searched and root mean square deviation (RMSD)-based top 20 naturally derived compounds were retrieved from the ZINC database [25] for docking analysis.

2.4. Molecular docking

The 2D structures of each chosen compound were manually designed using ChemDraw Professional 16.0. The 2D structures were then converted to 3D using Chem3D 16.0 [26], and the resulting files were saved as sdf. The Chem3D Gaussian interface was used to minimize energy. The SDF files of each ligand underwent ligand processing to assign the proper bond order [27].

The protein data library provides the 3D crystal structures of the target proteins for download. Protein preparation was carried out using Mgl tools 1.5.7 and Biovia Discovery Studio software in accordance with the usual procedure. Cofactors and water molecules were removed. After removing the previously joined ligands, AutoPrep was used to add polar hydrogens and Kollman’s charges to create the protein [28].

The “receptor cavity method” was employed to predict the binding sites of receptor proteins using BIOVIA Discovery Studio 2021. The SDB-Site module of the BIOVIA Discovery Studio program enabled the identification and description of protein structural binding sites. This process utilized the inhibitory properties of residues present in the binding sites with center-x, y, and z values [27].

The molecular docking was performed using AutoDock Vina 1.5.7 [29]. Following preparation, the receptor and ligand structures were saved in pdbqt format and transferred to the Vina folder. Vina’s AutoDock application was launched using the Command Prompt (CMD). The program command executed was “vina --config conf.txt --log log.txt” [30]. This method uses a grid-based model of protein-ligand potential interactions to calculate the binding affinity. Soft-core potentials, which are employed by Vina tools, are helpful in generating a range of random conformations of tiny organics and macromolecules inside the active area of the target protein. After docking the ligands to the proteins, the relative intensity of their interactions was assessed to identify possible treatment options. The test ligands were docked after the co-crystallized ligand was redocked to ascertain its binding affinity for the target protein [31]. Lastly, the docking results were displayed in terms of docking score, and top-identified compounds were visualized using PyMOL (version 2.4.0) for a 3D view [32] and BIOVIA Discovery Studio 2024 (v24.1.0.23298) for 2D interaction [33,34].

2.5. ADMET analysis

To predict the effectiveness and safety of top-identified compounds, their ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties were determined. A web server, SwissADME [35], was used to predict the ADME (absorption, distribution, metabolism, and excretion) properties of the compounds. The top hits for both targets were compared with reference drugs. Toxicity analysis was performed using StopTox [36], a web server. The Simplified Molecular Input Line Entry System (SMILES) string of the compounds was incorporated as input files to run predictions in both web servers.

2.5. Bioactivity prediction

For ligand-based target prediction, Swiss Target Prediction [37] was applied, which determined the most likely protein targets. The top-identified compounds for both targets, along with the reference drug, were subjected to bioactivity prediction.

2.6. Density function theory analysis

The Gaussian 09 program [38] was employed to perform density function theory (DFT) calculations for the reference drug and the top-identified compounds for each target, followed by visualizations of molecular electrostatic potentials (MEPs). For conductor-like polarizable continuum model (CPCM) optimization and frequency calculations for physiological phase, the B3LYP approach was selected. In addition, the HOMO–LUMO energy gap was optimized using Gaussian 09, with VESTA (version 3) applied to visualize their models for enhanced computational clarity [39,40].

2.7. Molecular dynamics simulations

Ligand-13 and reference drug in complex with Target-2 (CK2α) underwent molecular dynamics (MD) simulations with Amber17 to gain more knowledge about the stability of the chemical structure of the ligand when bound to the target protein. Initially, the process was facilitated by assigning partial atomic charges to Ligand-13 and the reference drug using the Antechamber module of Amber. Then, both compounds and the target protein (CK2α) were subjected to the Leap module, in which missing hydrogen atoms were added, the system was neutralized, and the complexes were solvated followed by the generation of necessary parameters and coordinate files for subsequent simulations. Then, for characterization of the dynamic properties of the system, the ff14SB force field was applied to the protein, while the generalized Amber force field (GAFF) was used for Ligand-13 and the reference drug. Two chloride ions were added to neutralize the protonated protein, and the complex was solvated in an octahedral box with TIP3P water, maintaining a margin distance of 10.0 Å. Finally, the necessary coordinate and parameter files were generated for the simulation, and the solvated complex was stored in the PDB format.

The complexes were reduced during the simulations to eliminate steric conflicts. Following minimization, the system was progressively heated from 0 to 300 K. To further stabilize the system, it underwent another equilibrium stage at 300 K. Equilibrium was established using Langevin dynamics with a force constant of 10 kcal/(mol Ų) and a collision frequency of 1 ps ⁻ ¹. After stabilization, the NPT ensemble was used to run 100 ns MD simulations at 300 K and 1 atm [41].

Several parameters, including RMSD, solvent-accessible surface area (SASA), and root mean square fluctuation (RMSF), were calculated at the end of the simulations using the CPPTRAJ module of Amber 17 [41] The radius of gyration was calculated for 100 ns using the method described by Arantes et al [42]. Principal component analysis and dynamic cross-correlation (DCCM) were also performed using the Bio3D package of R by writing a script specifically for these calculations [43–45].

2.8. Binding free energy MMGBSA calculations

Using MD shots, the binding free energy was computed to learn more about the energetic values and structure of the complexes. The MMPBSA/MMGBSA modules were used in conjunction with Amber17 to accomplish this [46]. These data were obtained by calculating the final 2 ns of data from 1000 snapshots of both complexes to determine the overall energies of Ligand-13 and the reference drug in complex with Target-2 (CK2α) using the following equation [47]:

$$[eqn]$$

The molecular mechanics energy (E_mm_) was calculated using the following equation:

$$[eqn]$$

The fragmentation of E_mm_ was analyzed as van der Waals (vdW) energies, non-bonded electrostatic energies (E_ele_), and solvation-free energy (G_sol_) encompassing both polar and non-polar:

$$[eqn]$$

3. Results

3.1. Target proteins

Fig 2 shows the 3D structures of the two target proteins after the removal of unnecessary associated structures such as ions and small molecules.

Selected target proteins of Plasmodium falciparum.A: cGMP-dependent protein kinase PKG. B: Protein kinase CK2 catalytic domain.

3.2. Pharmacophore-based virtual screening

The naturally-derived compounds from the PBVS were sorted based on their RMSD values, and the top 20 compounds with the lowest RMSD values were downloaded (Table 1) and visualized along with the reference drug in a ball-and-stick presentation using PyMOL (Fig 3).

Table 1: Top 20 Pharmit-screened naturally derived compounds arranged in increasing order of root mean square deviation (RMSD) values.

The 3D structures of the reference compound (U) and top 20 pharmacophore-based virtually screened naturally derived compounds (Ligand-1 to Ligand-20 presented chronologically as A to T).

3.3. Molecular docking

The docking results (Table 2) showed that the reference drug had a docking score of −9.054 kcal/mol and −7.683 kcal/mol with PKG and CK2, respectively. For PKG, the reference inhibitor showed the best docking score, followed by Ligand-9 with a score of −7.490 kcal/mol. For CK2, the best docking score was with Ligand-13, i.e., −11.468 kcal/mol. Thus, Ligand-9 and Ligand**-13** were identified as the lead compounds.

Table 2: Docking scores of target proteins (PKG and CK2) with the reference drug (PubChem CID: 164628676) and pharmacophore-based virtually screened naturally derived compounds generated by Maestro.

Ligand-9 and Ligand-13, the top compounds identified based on the docking score, were selected for their specific interactions with Target-1 (PKG) and Target-2 (CK2), respectively. The interactions were visualized on Discovery Studio and PyMol in comparison with the reference drug. All four selected docked complexes exhibited van der Waals, conventional hydrogen bond, and carbon-hydrogen bond interactions, while pi–cation, pi–sulfur, and salt bridges were observed selectively (Supplementary Table S1 in S1 File).

In the PKG docked complexes (Fig 4), the key residues involved in interactions were PHE-683, ASP-682, LEU-671, THR-622, VAL-621, LEU-620, PHE-616, ILE-602, THR-618, THR-593, LYS-570, LEU-569, LEU-557, VAL-555 and ILE-547. VAL-621 and LYS-570 contributed to strong hydrogen bonding, which enhanced structural rigidity of complex with compound 9. The pi–sigma interaction was found only in ILE-547 of the PKG-compound 9 complex and contributed to binding specificity and stabilization. ILE-681 and THR-618 stabilized the complex with reference drug by forming pi-sigma bonds. The pi–alkyl interactions that led to hydrophobic stabilization were observed at LEU-671, LEU-620 and VAL-555 in both complexes.

Visual representation of the site binding of Target-1 (PKG) to the reference drug and Ligand-9.A: 3D representation of PKG protein. B: PKG–Ligand-9 3D and 2D interaction with labeled interacting residues. C: PKG–reference drug 3D and 2D interaction with labeled interacting residues.

In the CK2 docked complexes (Fig 5), the main residues involved in interactions included VAL-57, ARG-51, GLY-50, GLY-52, ASN-165, MET-167, PHE-117, SER-179, LYS-72, ILE-178, ALA-70 and ILE-99. SER-55 formed strong H-bond with CK2. Interestingly, in the CK2–Ligand-13 docked complex, ARG-51 contributed to strong hydrogen bonding, while GLY-50 and ASN-165 participated in carbon-hydrogen bonding and SER-179 involved in pi-donor H-bond, enhancing the structural rigidity. PHE-117 and ILE-178 formed pi-sigma bond with ligand-13, while the pi-sigma interaction was found only in ILE-178 of CK2- reference drug complex contributing to the binding specificity and stabilization. In CK2-reference docked complexes, the pi–sulfur interactions were observed in MET-167. These pi–sulfur interactions resulted from the ability of sulfur atoms to engage with pi-electron clouds, which promoted hydrophobic and electrostatic complementarity between the receptor and ligand in the docked complex.

Visual representation of the site binding of Target-1 (CK2) to the reference drug and Ligand-13.A: 3D representation of CK2 protein. B: CK2–Ligand-13 3D and 2D interaction with labeled interacting residues. C: CK2–reference drug 3D and 2D interaction with labeled interacting residues.

3.4. ADMET analysis

This study analyzed the ADME of reference drugs and lead compounds (Ligand-9 and Ligand-13) and found an acceptable trend (Table 3). The compounds showed good gastrointestinal absorption and optimal oral bioavailability (Fig 6).

Table 3: ADME properties of the reference drug (A), Ligand-9 (B), and Ligand-13 (C).

Boiled egg diagram of the reference drug and Ligand-9 and Ligand-13, in which molecule-1 refers to the reference drug, molecule-2 refers to Ligand-9, and molecule-3 represents Ligand-13.The radar charts in the left corner show A: reference drug; B: Ligand-9; and C: Ligand-13.

Physiochemical properties indicated that the compounds did not violate major drug development guidelines, based on Lipinski’s rule of five. All the compounds were moderately soluble, with Ligand-13 showing the best solubility. The compounds also had moderate synthetic accessibility, with Ligand-9 being easier to synthesize than the reference drug. However, the compounds’ enzyme inhibition properties varied slightly, with Ligand-13 inhibiting the smallest number of enzymes, indicating lower drug interactions and efficient metabolism. Furthermore, the negative log Kp values of all compounds indicated their poor skin permeability, which can be an important consideration in their topical use. Moreover, Ghose, Veber, and Egan violations were also not observed for any compound, further strengthening their drug-likeness. There were also no PAINS alerts, which indicated the target specificity of the top-identified compounds and the reference drug.

In the toxicity analyses (Table 4), all three compounds showed eye irritation and corrosion. The reference drug and Ligand-13 showed oral toxicity, but acute inhalation toxicity was not found. Ligand-9 uniquely had a nontoxic oral toxicity prediction.

Table 4: StopTox toxicity parameters of the reference drug (A), Ligand-9 (B), and Ligand-13 (C).

3.5. Bioactivity

In the bioactivity prediction, kinase, enzyme, family A-G protein-coupled receptor, and cytosolic proteins were found to be common targets of all three compounds. Ligand-13 was found to be active against most biological targets (Fig 7C), indicating the diversity of its interactions. For the reference drug, the major targets were family A-G protein-coupled receptors, while kinase and enzyme proteins were secondary (Fig 7A). Ligand-9 (Fig 7B) was also predominantly active against family A-G protein-coupled receptors, followed by oxidoreductase and other cytosolic proteins.

Bioactivity of the reference drug (A), Ligand-9 (B), and Ligand-13 (C).

3.6. Density function theory calculations

The electronic properties and binding energies of the top-identified compounds, compared to those of the reference drug, revealed important details regarding their reactivity and stability (Table 5). All the compounds were predicted to be moderately reactive in terms of DFT-calculated parameters. The reference drug was found to be reactive and stable at the same time, owing to an energy gap of 4.59 eV and a high ionization potential and hardness. Its electrophilicity of 6.172 eV was comparatively moderate. Ligand-9 showed the highest chemical reactivity due to the smallest energy gap of 3.79 eV, the highest softness score of 0.527 eV, and moderate electrophilicity (6.038 eV). In contrast, its lower ionization potential (5.283 eV) marked it as the least stable of the three compounds. Ligand-13 was found to be intermediately reactive and stable due to an energy gap of 4.14 eV and 5.283 eV ionization potential, respectively. However, it had the highest electrophilicity (7.073 eV), electron affinity (1.759 eV), and electronegativity (3.827 eV) of all three compounds, indicating its increased reactivity due to increased capability of accepting electrons. Furthermore, it also had the highest dipole moment of 7.27 Debye, which indicated the strongest molecular polarity.

Table 5: Parameters of density function theory analysis of the reference drug (A), Ligand-9 (B), and Ligand-13 (C).

In addition, MEP surfaces were created to predict the regions of the selected compounds that might participate in electrophilic or nucleophilic reactions. The reference drug had an even potential distribution, while Ligand-9 and Ligand-13 had strong potential sites as indicated by the color-coded MEP in Fig 8.

Visual representation of the molecular electrostatic potential and HOMO–LUMO energy gap of the reference drug (A), Ligand-9 (B), and Ligand-13 (C).

Ligand-13 was just above the classification threshold in the StopTox model, with an estimated 55% chance of acute oral toxicity. Given that substances with a high resemblance to poisonous references usually score ≥70–80%, this suggests a borderline or moderate risk rather than a substantial toxic liability. Ligand-9 performed worse than Ligand-13 in crucial measures, such as the CYP enzyme inhibition profile, DFT reactivity descriptors, and docking affinity, which are crucial for drug development, even though Ligand-9 was projected to be non-toxic (51%). For MD simulation, Ligand-13 was prioritized for additional evaluation in the MD simulation.

3.7. Molecular dynamics simulations

Ligand-13 in complex with Target**-2 (CK2)** was narrowed down for MD simulation due to stronger binding affinity and maximum stability with most biologically active targets. For computational ease, only the actively interacting chain of the target protein (CK2α) was included in the simulations. The RMSD graphs (Fig 9) revealed overall stable binding of the docked complexes as indicated by a plateauing trend in the graphs. In both complexes, the target protein, their binding site, and the bound compounds, were all in close proximity. The deviations of the docked complex with the reference drug were within an acceptable range, i.e., 3.5 Å, and were prominently different around 16 ns and 54 ns. On the other hand, the docked complex with Ligand-13 showed that all the deviations were notably different around 16 ns and 54 ns and ranged within 2.5 Å.

Root mean square deviation (RMSD) of the docked complex of CK2 with reference drug (A) and with Ligand-13 (C), with the highlighted superimposed structures B and D, respectively.The RMSD plot of the protein–ligand complex displayed two notable fluctuations around 16 ns and 54 ns, indicating transitions between distinct conformational states. Based on the RMSD trajectory, the system remained stable during 0–16 ns (Basin I), shifted to a new conformation from 16–54 ns (Basin II), and then returned toward the initial stable state after 54 ns (Basin I).

The RMSF (Fig 10A) showed a similar trend, in which slight fluctuations were observed around residue numbers 100 and 260 of the protein, which were attributed to the terminal region of the protein structure. The radius of gyration of the docked complexes was within 21.2 Å, indicating that they remained compact throughout the simulation (Fig 10B). The SASA varied between 400 and 450 Å, suggesting that the protein structures underwent minute conformational changes when exposed to solvents (Fig 10C).

A, B, and C plots of the docked complex of CK2 with Ligand-13 (blue) and the reference drug (purple) show the root mean square fluctuation (RMSF), radius of gyration (Rg), and solvent-accessible surface area (SASA), respectively.

During simulations with the reference medication, principal component analysis revealed compact clustering of the protein confirmations, indicating the stability of their docked complexes (Fig 11A, 11B). However, the highly uneven distribution of the CK2-Ligand-13 complex suggests that it becomes more flexible upon binding.

Comparative representation of principal component analysis (PCA), free energy landscape, and dynamic cross-correlation of the docked complex of CK2 with the reference drug (A, C, E) and with Ligand-13 (B, D, F).

The free energy landscapes (FEL) of both docked complexes (Fig 11C, 11D) showed deep blue basins at various points, indicating energetically stable conformational states during ligand binding. In accordance with the conformational variability of the complexes, the large red areas were clearly indicative of higher-energy conformations. Interestingly, the FEL showed two different minima, suggesting that both complexes had two stable conformational states during the simulations. The shallower blue region, observed between 16 and 54 ns, represents a momentarily less stable state before the system returns to its original basin, whereas the deep blue region corresponds to Basin I, which is the most stable conformation observed at 0–16 and after 54 ns. The complex underwent reversible structural transformations prior to settling in an energetically advantageous state, as indicated by the consistency between the RMSD and FEL. This conclusion is further supported by the quantitative examination of basin populations and average RMSD/FEL values (Supplementary Tables S2 and S3 in S1 File).

Finally, the correlated and anticorrelated movements of the protein in both docked complexes were revealed by a consistent pattern inter-residue DCCM analysis (Fig 11E, 11F), with positive and negative correlations indicated by red and blue colors, respectively. It was noteworthy that for the residues 0–100 (N-lobe of CK2 chain A), the negative correlations were dominant in Ligand-13’s complex. In contrast, Ligand-13 exhibited increased positive correlations from residues 100–300 (C-lobe of CK2α).

During simulations, the directional movements of the docked complexes in their minimum energy confirmations were also visualized (Fig 12). These observed movements were compared to the protein structure (CK2α). The shorter N-lobe of the protein structure, which also contained its active binding site and P-loop, was observed to constantly display movements throughout the simulation. On the other hand, the C-lobe also displayed movements, but in the opposite direction to the N-lobe. These movements were significant for the effective docking of compounds.

Comparative visualization of the directional movements of protein during the docking of CK2 with the reference drug (A, B) and with Ligand-13 (C, D) observed at different intervals.The detailed reference structure of the protein is presented in the center.

3.8. Free binding energy MMGBSA calculations

In the comparative analysis of the docked complexes of the reference drug and Ligand-13, both complexes’ ΔE and ΔG showed similar results, except for a large difference in ΔG values in the gas and solvent phase shown in Fig 13 (Supplementary Table S4 in S1 File). The docked complex of the reference drug exhibited more negative ΔGgas signaling toward its strong gas-phase interactions, which were also supported by van der Waals and electrostatic interaction energy values. On the other hand, the docked complex of Ligand-13 displayed better binding in solvated states, as predicted by its ΔE and ΔG values*.*

Comparative analysis of ΔE and ΔG values of the reference drug (purple) and Ligand-13 (blue).

4. Discussion

In this study, we explored the potential of naturally derived compounds as novel antimalarial agents targeting P. falciparum’s protein kinases, PKG and CK2, which are essential to various stages of the parasite’s life cycle [14,21]. Given the rise in resistance to antimalarial treatments, natural compounds and their derivatives are important in drug discovery. Natural materials have traditionally offered a diverse chemical pool for the development of antimalarial drugs [48]. Plants are the source of several effective antimalarial medications, such as quinine [49] and artemisinin, which form the basis of contemporary therapy. However, there has been a rise in pathogen resistance to these medications [50]. This demonstrates how the natural world may offer the next wave of treatments for malaria that can successfully combat drug resistance.

With docking scores of −7.490 kcal/mol and −11.468 kcal/mol of Ligand-9 and Ligand-13 respectively, were identified as strong inhibitors of the target proteins (PKG and CK2) in this study. The strong binding affinities of these substances, particularly Ligand-13, indicate their considerable capacity to disrupt the life cycle of P. falciparum. The significance of PKG and CK2 in parasite pathogenicity was confirmed by their suppression. This study expands on earlier research on the inhibition of CK2 and PKG, which have been shown to be effective therapeutic targets for the treatment of malaria [12,51].

Based on the outcomes of specific interactions, including pi–sulfur interactions, hydrogen bonding, van der Waals forces, and special salt bridges, these compounds can effectively inhibit target proteins by forming stable complexes with them. Other studies have also confirmed the inhibition of the resistance of particular P. falciparum through compounds with similar docking scores as observed in this study [21]. These results underscore the significance of these naturally derived products as potential candidates for the development of antimalarial drugs.

The ADME properties of Ligand-9 and Ligand-13 were found to be acceptable enough to consider them good antimalarial drug candidates. Gastrointestinal absorption, important for the efficacy of the drug in the human body, was found to be high for both these compounds. However, they were not predicted to cross the blood–brain barrier (BBB), which should be addressed specifically in relation to the disease and the mechanism of action of drugs. Generally, P. falciparum refrains from targeting brain cells. However, in cases of faulty diagnosis and improper medication, this parasite may invade brain cells, giving rise to a severe progressive form of the infection called cerebral malaria [52]. The naturally derived compounds discussed in this study are proposed to be effective for the widespread form of P. falciparum infections. However, they can be used against the more severe form of the disease if they are used in conjugation with other BBB-crossing antimalarial compounds [53]. The predicted bioactivity of Ligand-9 and Ligand-13 further reinforces the use of these compounds as antimalarials. Ligand-13, in particular, was found to be effective against multiple biological targets and might be highly effective in targeting the complexity of malarial infection [54].

Concerns regarding the safety of Ligand-13 and the reference compound as antimalarials are raised by their expected oral toxicity. However, toxicity is dose-dependent and can be reduced by changing the concentration or dosage of the drug [55]. Oral toxicity can be readily reduced by efficient dose adjustment of these drugs. Additionally, the distribution of medications can be modified to ensure their safe and efficient use [56]. Ligand-13’s intermediate oral toxicity prediction of 55% does not hinder its advancement in therapeutic development. Given its superior pharmacological profile and the hope that future formulations or dosing schemes would successfully reduce its mild oral toxicity signal, Ligand-13’s priority is still scientifically warranted.

For optimal drug designing, DFT-calculated parameters are essential [57]. For this purpose, the electronic properties and reactivity of Ligand-9 and Ligand-13, in comparison with the reference drug, were determined. This analysis showed that Ligand**-9** with the lowest energy gap (3.79 eV) was the most reactive but it was less stable due to its lower ionization potential, as indicated by similar results in other studies [58]. However, Ligand-13 showed optimum drug properties with an intermediate HOMO–LUMO energy gap (4.14 eV), maximum electrophilicity, and a dipole moment of 7.27 Debye that provides suitable interaction ability.

As docking is a static representation of a molecule’s binding position in the protein’s active site, MD simulations often integrate Newton’s classical equation of motion to calculate atom motions over time [59]. Molecular dynamics simulations were used to predict the ligand-binding status in a physiological environment. MD simulations of the reference drug and Ligand-13 in complex with CK2 revealed the stable binding of protein with both these compounds. Our 100 ns MD simulation findings are consistent with earlier research that was carried out over longer timeframes (200 ns), which also demonstrated a robust and consistent interaction between the amino acids and the inhibitors of the simulation’s duration [60,61]. The conformational changes observed in the protein structure over time were similar for both docked complexes, which are predictive of similar inhibition of the target protein by Ligand-13. The movements observed in the N-loop of the protein signal toward flexible accommodation of Ligand-13 in the active site are crucial for its inhibitory functions. This region is particularly significant for the functioning of CK2; its inhibition can result in the elimination of P. falciparum [24]. Our study reaffirms that naturally derived products, particularly Ligand-13, hold promise for overcoming the resistance to current therapies.

In contrast to earlier in silico screenings of P. falciparum kinases, this study combined multi-target parallel docking against PKG and CK2 with large-scale pharmacophore-based screening of a ZINC natural subset. This was followed by DFT reactivity assessment, MD-based stability analysis, and MM/GBSA binding energy profiling. This integrated workflow prioritizes Ligand-13 as a lead for additional biological assessment by reducing the number of candidates based on projected reactivity, safety liability, and binding affinity.

Although in silico tools highlighted useful naturally derived antimalarial compounds, these findings need to be validated via in vitro and in vivo experiments. In addition, more drug targets of P. falciparum can be used to combat the multi-targeted drug resistance of this deadly pathogen.

5. Conclusion

By assessing the activity of naturally occurring chemicals, our in silico method proposed Ligand-13 as a possible antimalarial medication candidate against probable P. falciparum targets, PKG and CK2. Comparing Ligand-13 to the other studied ligands, Ligand-13 showed the most promising combination of robust electronic characteristics, good CYP enzyme interaction profile, and greater binding affinity. Together, our results provide valuable new insights into the potential of Ligand-13 as a lead scaffold. Ligand-13 has the potential to be a lead compound, even though in silico models indicate a moderate probability of acute oral toxicity (55%) for this compound. Instead, it emphasizes the importance of continuing to optimize the medication development process. Future research should focus on dose-dependent toxicity assessment, alternative administration methods, and structural refinement to reduce toxicity risk. Confirming the safety profile of the compound would also require experimental validation using in vitro cytotoxicity and metabolic stability assays, in vivo acute oral range-finding, and hERG liability investigations.

Supporting information

S1 FileSupplementary data related to docking and molecular dynamics analyses.This file contains Supplementary Tables S1–S4, including the summary of interacting residues and interaction types in docked complexes, quantitative RMSD and free energy landscape (FEL) analyses from 100 ns molecular dynamics simulations of Ligand-13 and the reference drug, and energy difference values of docked complexes with target 2.(DOCX)