Physics-based simulations applied to macromolecules: Plasmodium falciparum Dihydroorotate Dehydrogenase (6 VTY) in complexation with a series of pyrrole-derived inhibitors

 Corresponding author: Niaré Adama Copyright © 2025 Author(s) retain the copyright of this article. This article is published under the terms of the Creative Commons Attribution License 4.0. Physics-based simulations applied to macromolecules: Plasmodium falciparum Dihydroorotate Dehydrogenase (6 VTY) in complexation with a series of pyrrolederived inhibitors Rika Justin Kouadja 1, Adama Niaré 1, *, Matthias Logbo Moussé 1, Ibrahim Bamba 1, Mamadou Guy-Richard Koné 2 and Eugène Megnassan 1, 3, 4, 5 1 Laboratory of Fundamental and Applied Physics, UFR SFA, Nangui Abrogoua University, BP 801 Abidjan 02, Côte d’Ivoire. 2 Laboratory of Thermodynamics and Physico-chemistry of the Environment, UFR SFA, Nangui University 02 BP 801 Abidjan 02, Côte-d’Ivoire. 3 Laboratory of Materials Science, Environment and Solar Energy, UFR SSMT, Félix Houphouët-Boigny University, BP 582 Abidjan 22, Côte d’Ivoire. 4 Laboratory of Structural and Theoretical Organic Chemistry, Félix Houphouët-Boigny University, BP 582 Abidjan 22, Côte d’Ivoire. 5 ICTP-UNESCO, QLS, Strada Costiera, 11, I-34151, Trieste, Italy. World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 Publication history: Received on 08 August 2025; revised on 14 September 2025; accepted on 16 September 2025 Article DOI: https://doi.org/10.30574/wjarr.2025.27.3.3231 Abstract Despite the efforts and resources devoted to combating malaria, as well as the knowledge acquired about the different species of Plasmodium, of which Plasmodium falciparum is the most common in humans, malaria remains the leading endemic parasitic disease in the world. Malaria alone is responsible for approximately 597,000 deaths. It mainly affects pregnant women and children under the age of 5. The emergence of drug-resistant strains of the parasite poses a major challenge in the treatment of malaria. The enzyme dihydroorotate dehydrogenase (DHODH, pdb code: 6VTY) plays a crucial role in the treatment of malaria through the pyrimidine biosynthesis pathway, which is essential for the growth and reproduction of the Plasmodium falciparum (pf) parasite. Physics-based simulations applied to macromolecules are a major asset in the development of new drugs capable of overcoming resistance issues, improving selectivity, and increasing clinical efficacy. Using structure-based molecular design, we have developed Quantitative Structure-Activity Relationship (QSAR) models of 3D complexation from pyrrole (Pyl) molecule derivatives. A linear correlation has been established between Gibbs free energies (GFE: ∆∆Gcom ) and the observed enzyme inhibition constants (IC50 exp) for each complex formed: pIC50 exp = - 0.1207×∆∆Gcom+ 7.8526; R2 = 0.95. The effect of the biological (aqueous) environment was taken into account as a component of ∆∆Gcom through the free energy of molecular electrostatic solvation (∆Gsol ) solution of the Poisson-Boltzmann equation. The predictive power of the RQSA model was validated by the generation of 3D-RQSA Pharmacophores (PH4): pKi exp = 1,1585×pKi pre − 1.1448; R2 = 0.93. The combination of 3D complexation and 3D-RQSA pharmacophore methodologies has provided key insights into the orientation and configuration of pyrrole (Pyl) molecules within the active site of dihydroorotate dehydrogenase (DHODH) for effective malaria control. Keywords: QSAR Model; Complexation Model; Pharmacophore Model; Molecular Modeling; Dihydroorotate Dehydrogenase; Pyrrole Family World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 962 1. Introduction Physics-based simulations have significantly increased the efficiency of computer-aided drug design (CADD). Molecular modeling (MM) and computer-aided combinatorial chemistry (CACC) have become new techniques for understanding many chemical, biological, and physical phenomena. These tools are increasingly useful for the pharmaceutical industry and enable the exploitation of natural substances [1]. The behavior of molecules can be simulated statically or dynamically [2]. This simulation includes, in particular, the behavior and flexibility of molecules during inhibition. The reaction mechanisms and interactions involved in enzyme catalysis are also included throughout this simulation [2]. Furthermore, physics-based simulation highlights complexation and the prediction of properties and functionalities based on structural analogies [2]. These techniques are proving to be essential in the development of new drugs to combat certain endemic diseases such as malaria, tuberculosis, HIV/AIDS, dengue fever, etc. These techniques are proving to be more than necessary in the development of new drugs to combat certain endemic diseases such as malaria, tuberculosis, HIV/AIDS, dengue fever, etc. In this paper, we focus on malaria due to the alarming statistics from the World Health Organization (WHO) [3]. Approximately 246 million people in 109 countries are exposed to malaria each year. Malaria remains the most significant parasitic disease, with more than 597,000 deaths, 77% of which are children under the age of 5 [3]. The situation has become a cause for concern because in most regions of the world, Plasmodium resistance has rendered traditional antimalarial drugs (quinine and chloroquine) ineffective, jeopardizing the fight against malaria [4]. Given this urgent need, the discovery and development of new therapeutic targets and new molecules are therefore required in order to expand the anti-malarial arsenal. Dihydroorotate dehydrogenase (DHODH), located in the mitochondria of the Plasmodium falciparum parasite and essential for its replication and development in red blood cells, is therefore an essential target for the treatment of malaria [5]. Pyrrole-based compounds are a diverse class of chemical compounds that target and block the activity of the enzyme dihydroorotate dehydrogenase (DHODH) [5]. In short, these pyrrole-based inhibitors play an essential role in the inhibition of DHODH [5]. The objective of this computer simulation work is to understand how pyrrole-based molecules interact within the active site of DHODH in order to evaluate the different binding forces that result. To achieve our objective, we developed a Quantitative Structure-Activity Relationship (QSAR) model based on DHODH-Pylx complexes obtained by in situ modifications of the crystallographic structure of the DHODH-Pyl1 complex [6, 7]. This model allowed us to establish a linear relationship between the relative change in Gibbs free energy of the formation of the DHODH-Pylx complex and inhibitory activities. This enabled us to identify the active conformation of Pylx at the active site of DHODH through complexation using Poisson-Boltzmann Molecular Mechanics (PB-MM). Using this active conformation, we created a 3D RQSA pharmacophore for DHODH inhibition. 2. Material and methods 2.1. Training and validation sets In this work, the molecules used are from the pyrrole family. They have very interesting inhibitory properties in studies of dihydroorotate dehydrogenase (DHODH) enzyme inhibition. These pyrrole derivative molecules were synthesized and biologically evaluated by Sreekanth Kokkonda et al [7]. Their inhibitory potencies cover a wide range of inhibitory concentrations between 0.014 μM and 1.6 μM (Table 2), enabling the design of an RQSA model. Out of a total of 16 molecules, this pyrrole (Pyl) derivative inhibitors were divided into two groups comprising 12 for the test set and 4 for the validation set using the “Generation training and test set” protocol integrated into the Discovery Studio software [8]. 2.2. Model construction Three-dimensional (3D) models of free inhibitors (I), free DHODH enzyme (E), and enzyme-inhibitor complexes (E: I) were created based on the high-resolution crystal structure (1.18 Å) of a reference complex with the inhibitor Pyl 1 (PDB code: 6VTY) [7]. These models were generated using the graphical interfaces of the Insight-II [9] and Discovery Studio 2.5 [8] molecular modeling tools. 2.3. Molecular Mechanics The modeling of DHODH and Pylx ligand complexes was performed using molecular mechanics, applying the CFF force field [10]. This process involves optimizing the geometry of ligand-enzyme interactions by minimizing energy, allowing for a more accurate representation of the binding and structural characteristics of the complexes. World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 963 2.4. Conformational search The conformations of free inhibitors were derived from their bound conformations within protein-ligand (PL) complexes. This was achieved by implementing a stepwise relaxation process, allowing the structures to gradually adjust until they reached the nearest local energy minimum [11]. 2.5. Solvation of Gibbs free energies Ligand-receptor interactions occur in a solvent environment, which affects the binding process through hydrogen bonding and solvation effects. To account for these influences, the electrostatic component of Gibbs free energy was calculated by solving the nonlinear Poisson-Boltzmann equation [12], incorporating ionic strength. This equation can be written as: ∇2V+ρ ε= 0 (1) This calculation was performed using the Delphi module implemented in Discovery Studio 2.5 [8]. 2.6. Calculation of binding affinity and QSAR model For each molecule in the test set, the relative thermodynamic value (ΔΔGcom) was calculated. This calculation was performed using the MM-PB molecular complexation method. Since the Pyl1 ligand was the most active (IC50 exp = 0.014 μM), it was chosen as the reference ligand. The models consisted of first correlating the variation in experimental biological activity (pIC50 exp) with that of the enthalpy (ΔΔHMM) of Molecular Mechanics (gas model) and then doing the same between the variation in experimental biological activity (pIC50 exp)and the relative variation in free enthalpy of complexation (ΔΔGcom) (solvent model). pIC50 exp was obtained from the following equation: 𝑝IC50 exp = −log10 (𝐼𝐶50 106 ⁄) (2) To validate the MM-PB molecular complexation 3D-RQSA model, statistical indicators of the reliability of linear regressions were calculated. 2.7. Pharmacophore generation The 3D-QSAR pharmacophore generation protocol (PH4) in Discovery Studio [8], based on the Catalyst HypoGen algorithm [13, 14], was used to develop the PH4 model for DHODH inhibition. 3. Results and discussion 3.1. Training and validation sets Throughout this work, 12 Pylx inhibitors were used for the test set and 4 Pylx inhibitors were used for the validation set, as shown in Table 1. The experimentally determined inhibitory activities are known and come from a single laboratory [7]. The entire series of inhibitors was obtained by variations at two positions, Ar and R, on the pyrrole nucleus skeleton, as shown in Table 2. The maximum experimental inhibitory concentrations (0.014 μM ≤ IC50 exp≤ 1.6 μM) [7] cover a sufficiently wide range of concentrations to allow the construction of a reliable QSAR model (Table 2). The ratio between the sizes of molecules in the test set and validation set remains a critical factor for correct classification, but it is limited by the number of sets of homologous compounds available in the literature [7, 15]. World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 964 Table 1 Training set (Pyl1-12) and validation set (Pyl13-16) of PDHODH inhibitors [7] used in the preparation of QSAR models of inhibitor binding. The Ar and R groups are numbered in the first part of the Table as #R ≡ group index. #R 1 2 3 4 R group #R 5 6 7 8 R group #R 9 10 11 12 R group #R 13 14 World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 965 Table 2. inhibitory activities (𝐼𝐶50 𝑒𝑥𝑝) of molecules used for the training set (Pyl1-12) and validation set (Pyl13-16) Training Set Pyl 1 Pyl 2 Pyl 3 Pyl 4 #Ar-#R 8-1 8-2 9-1 13-5 𝐼𝐶50 𝑒𝑥𝑝 (μM) 0.014 0.019 0.021 0.022 Training Set Pyl 5 Pyl 6 Pyl 7 Pyl 8 #Ar-#R 13-3 6 10-1 7-1 𝐼𝐶50 𝑒𝑥𝑝 (μM) 0.031 13-6 0.040 0.045 Training Set Pyl 9 Pyl 10 Pyl 11 Pyl 12 #Ar-#R 12-3 11-1 11-6 2 𝐼𝐶50 𝑒𝑥𝑝 (μM) 0.059 0.036 0.220 12-4 Validation Set Pyl 13 Pyl 14 Pyl 15 Pyl 16 #Ar-#R 12-5 14-1 13-2 13-3 𝐼𝐶50 𝑒𝑥𝑝 (μM) 0.023 0.036 0.036 1.600 3.2. QSAR Model 3.2.1. One Descriptor QSAR Models The complexation of the 12 inhibitors in the test set and the 4 inhibitors in the validation set was achieved by in situ modification of the crystal structure of complex 1, in this case the first inhibitor in the test set series (DHODH: Pyl1), following the procedure detailed in the Methods section. The relative Gibbs free energy (GFE) of complex formation (DHODH: Pyl) was calculated for each of the 16 optimized (enzyme-inhibitor) complexes. Table 3 presents the calculated energy values and their components. To establish a quantitative structure-activity relationship (QSAR), the experimental inhibitory activities (IC50 exp) of the Pylx compounds were converted to (𝑝IC50 exp) using the using equation 2. The linear relationship derived from this RQSA model between pIC50 exp and ΔΔGcom is highlighted in the following equation: pIC50 exp = - 0.1207×∆∆Gcom+ 7.8526 (3) The resulting strong correlation confirmed the binding modes of Pylx within the DHODH site, leading to the definition of the DHODH inhibition pharmacophore. In order to better understand the binding affinity of Pylx to DHODH, the gasphase complexation enthalpy (ΔΔHMM) was calculated and analyzed. The statistical significance of this linear relationship (Table 4, Equation A) highlighted the crucial role of inhibitor-enzyme interactions, particularly when solvent effects and entropy loss upon binding to DHODH were neglected. This correlation explained approximately 90% of the variation in values, highlighting the contribution of enthalpy to ligand binding affinity. In addition, a more comprehensive descriptor, the Gibbs free energy (GFE) of the formation of the Pyl1 complex, integrating all key components ΔΔHMM, ∆∆TSvib and ∆∆Gsol, was evaluated (Table 4, equation B). The relatively high correlation coefficient (R² = 0.95) demonstrated that structural information derived from 3D models of the complexes (DHODH: Pylx) could reliably predict how pyrrole inhibitors interact within the active site of DHODH using QSAR model B, as shown in Table 3. These values (Table 4) indicate that there is a strong relationship between the bonding model and the experimental model of the pyrrole series. World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 966 Table 3 Gibbs free energy (binding affinity) and its components for the training set of DHODH inhibitors Pyl1-12 and validation set inhibitors Pyl13-16 [7] Training set a MW b ∆∆HMM c ∆∆Gsol d ∆∆TSvib e ∆∆Gcom f 𝑰𝑪𝟓𝟎 𝒆𝒙𝒑 g [g.mol-1] [kcal.mol-1] [kcal.mol-1] [kcal.mol-1] [kcal.mol-1] [μM] Pyl 1 311.299 0.00 0.00 0.00 0.00 0.014 Pyl 2 325.326 1.10 -0.47 -1.38 2.01 0.019 Pyl 3 324.341 2.88 -1.68 0.27 0.93 0.031 Pyl 4 351.366 1.23 0.85 -1.38 3.45 0.036 Pyl 5 285.381 1.24 1.38 0.06 2.55 0.04 Pyl 6 33.34 3.09 0.67 -1.18 4.95 0.022 Pyl 7 287.354 3.32 -0.75 -0.06 2.63 0.045 Pyl 8 326.314 2.28 1.86 0.01 4.13 0.06 Pyl 9 287.354 2.26 2.42 -0.97 5.65 0.021 Pyl 10 338.367 1.34 0.06 -1.99 3.39 0.059 Pyl 11 351.366 7.31 -0.25 -1.81 8.86 0.22 Pyl 12 352.394 16.03 3.24 1.74 1.53 1.8 Validation set a MW b ∆∆HMM c ∆∆Gsol d ∆∆TSvib e ∆∆Gcom f 𝑝𝐶50 𝑝𝑟𝑒/𝑝𝐶50 𝑒𝑥𝑝 h [g.mol-1] [kcal.mol-1] [kcal.mol-1] [kcal.mol-1] [kcal.mol-1] Pyl 13 322.325 0.13 2.55 -3.02 5.71 0.94 Pyl 14 329.289 2.44 -0.16 -2.51 2.63 0.99 Pyl 15 325.329 3.11 1.64 -1.49 6.24 0.95 Pyl 16 339.355 5.13 1.29 -1.61 8.86 0.91 a for the chemical structures of the training set of inhibitors see Table 1; b Mw is the molar mass of inhibitors; c ∆∆HMM is the relative enthalpic contribution to the GFE change related to E-I complex formation derived by MM;∆∆HMM ≈[EMM{E-Ix} − EMM{Ix}] − [EMM{E-Iref} − EMM{Iref}], Iref is the reference inhibitor Pyl1; d∆∆Gsol is the relative solvent effect contribution to the GFE change of E-I complex formation: ∆∆Gsol = [Gsol{E-Ix} − Gsol{Ix}] − [Gsol{E-Iref} − Gsol{Iref}]; e ∆∆TSvib is the relative entropic contribution of inhibitor to the GFE of E-Ix complex formation: ∆∆TSvib = [TSvib{Ix}E − TSvib{Ix}] − [TSvib{Iref}E − TSvib{Iref}]; f∆∆Gcom is the overall relative GFE change of E-Ix complex formation: ∆∆Gcom ≈∆∆HMM + ∆∆Gsol − ∆∆TSvib; g IC50 exp is the experimental half-maximal inhibition concentration of PfA.M1obtained from ref. [Error! Bookmark not defined.]; h ratio of predicted and experimental half-maximal inhibition concentrations pC50 pre/pC50 exp (pIC50 pre = −log10 IC50 pre ) was predicted from computed ∆∆Gcom using the regression equation for DHODH shown in Table 4, B. World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 967 Table 4 Analysis of computed binding affinities ΔΔGcom, its enthalpic component ΔΔHMM, and experimental half-maximal inhibitory concentrations 𝐩𝐈𝐂𝟓𝟎 𝐩𝐫𝐞 = −log10 𝐈𝐂𝟓𝟎 𝐩𝐫𝐞 of Pylx towards DHODH [7] Statistical Data of Linear Regression A B pIC50 exp = - 0.1263 × ΔΔHMM + 7.7312 (A) pIC50 exp = - 0.1207 × ∆∆Gcom + 7.8526 (B) Number of compounds n 12 12 Squared correlation coefficient of regression R2 0.90 0.95 LOO cross-validated squared correlation coefficient R2xv 0.89 0.94 Standard error of regression  0.19 0.14 Statistical significance of regression, Fisher F-test 94.20 184.04 Level of statistical significance  >95% >95% Range of activities 𝐼𝐶50 𝑒𝑥𝑝 [µM] 0.014 – 1.6 0.014 – 1.6 The statistical data confirmed validity of the correlation Equations (A) and (B) plotted on Figure 1. The ratio pC50 pre/pC50 exp ≈1 (the pIC50 pre values were estimated using correlation Equation B, Table 4) calculated for the validation set Pyl13-16 documents the substantial predictive power of the complexation QSAR model from Table 2. A B Figure 1 (A) plot of correlation equation between 𝒑𝑰𝑪𝟓𝟎 𝒆𝒙𝒑 and relative enthalpic contribution to the GFE (∆∆HMM [kcal.mol-1]). (B) similar plot for relative complexation Gibbs free energies of the DHODH-Pyl complex formation ∆∆Gcom [kcal.mol-1] of the training set [7]. The validation set data points are shown in red color World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 968 3.3. Binding Mode of Pyls Beyond the robustness of the QSAR model, analysis of interactions between Pylx and DHODH active site residues within a radius of 10 Å revealed the key interactions responsible for the affinity of pyrrole inhibitors (Pylx) with DHODH. These interactions were highlighted through the Pyl 1 inhibitor, the most active inhibitor in the test set, as shown in Figure 2. Van der Waals interactions, low-intensity electrical interactions, were observed between Pyl1 and the residues Tyr168, Phe 171, Leu 172, Gly 181, Arg 265, and Gly 535 of DHODH. A T-shaped Pi-Pi stacking interaction, which is attractive and non-covalent, was observed with the Phe168 residue of DHODH. In addition, a Pi-Sulfur interaction was observed between the Cys 175 residue of DHODH and the phenyl ring of Pyl 1. Pi-alkyl interactions were observed between residues His 185, Phe 227, Ile 263, and Tyr 528 of DHODH with Pyl1. Hydrogen interactions were observed in small quantities in almost all complexes. These interactions are important for protein folding, conformational changes, and protein/ligand recognition. On the ligand side, they affect the physicochemical properties of molecules, such as solubility and membrane permeability, which are crucial elements in drug development. A C Figure 2 (A) 2D interaction diagram of the most active inhibitor (Pyl 1) in the training set within the active site of DHODH, (B) Schematic 3D interaction of Pyl 1 at the active site of the DHODH enzyme, (C) Hydrophobic surface of the active site of DHODH. Surface coloring legend: red ≡ hydrophobic, blue ≡ hydrophilic and white ≡ intermediate B World Journal of Advanced Research and Reviews, 2025, 27(03), 961-972 969 3.4. 3D-QSAR Pharmacophore Model 3.4.1. PfA-M1 Active Site Pharmacophore The interaction generation protocol in Discovery Studio molecular modeling program [8] provides the pharmacophore features of the active site of a protein. DHODH predominantly displays aliphatic and hydrophobic aromatic features at the active site. 3.4.2. Generation and Validation of 3D-QSAR Pharmacophore A 3D-QSAR pharmacophore model for DHODH inhibition was constructed based on the active conformations of the 12 Pylx inhibitors in the training set (TS) and validated by the 4 inhibitors in the validation set (VS). The 16 pyrrole-derived Pylx inhibitors cover a wide range of experimental activity (0.014 µM - 1.6 µM) (Table 2). The 3D-QSAR pharmacophore generation process is divided into three main phases: the constructive phase (i); the subtractive phase (ii); and finally, the optimization phase (iii). The hypotheses were scored based on errors in the activities estimated from regression and complexity using a simulated annealing approach. At the end of the optimization, the 10 highest-scoring unique pharmacophore hypotheses were retained, all displaying three characteristics: cost values, correlation coefficients, and mean squared deviation values (RMSD). The pharmacophore characteristics and maximum fit values of the 10 highestranked hypotheses (Hypo1-Hypo10) are shown in Table 5. They were selected based on significant parameters, such as a high correlation coefficient, low total cost, and low RMSD. Their selection was based on significant statistical parameters, including a high correlation coefficient, low total cost, and low RMSD. The cost difference (Δ = 184.5171) was calculated as follows: zero cost (216.723) – fixed cost (32.2059), indicating a high probability (> 90%) that the model represents a true correlation. Standard parameters such as RMSD vary between 1.72 and 4.46, and the correlation coefficient squared (R2) ranges from 0.96 to 0.63. The first hypothesis of PH4 was selected for further analysis. Its evaluation consisted of mapping the most active ligand, Pyl1, from the set of 12 inhibitors in the test set (Figure 3, B), showing its active conformation with the Hypo1 pharmacophore for DHODH inhibition. The statistical data for all hypotheses (costs, RMSD, R2) are listed in Table 5. The linear correlation equation of the 3D-QSAR pharmacophore model explaining the variation in experimental biological activity as a function of predicted activity is given by the following linear relationship: pIC50 exp = 1.1585×pIC50 pred −1.1488 (4) This relationship is illustrated in Figure 4E with its internal validation reliability statistical parameters such as n = 12, R2= 0.93, 𝑅𝑋𝑉 2= 0.91, F-test = 123.43, α = 0.16 < 1, with a confidence level of 95%. These results indicate that the 3DQSAR pharmacophore model is highly reliable and can be used to search for new competitive analogues for DHODH inhibition. Table 5 Parameters of 10 generated PH4 pharmacophoric hypotheses for DHODH inhibitor after Cat-Scramble validation procedure (49 scrambled runs for each hypothesis at the selected level of confidence of 98%) Hypothesis RMSD a R2 b Total Costs c Costs Difference d Closest Random e Hypo1 1.72 0.96 63.10 153.61 45.79 Hypo2 2.21 0.92 64.09 152.62 50.47 Hypo3 2.79 0.87 82.42 134.29 52.706 Hypo4 3.73 0.76 121.82 94.90 53.47 Hypo5 3.78 0.75 124.22 92.50 54.77 Hypo6 3.81 0.75 130.66 86.06 57.87 Hypo7 4.00 0.71 131.23 85.48 58.38 Hypo8 4.10 0.70 136.05 80.67 62.63 Hypo9 4.13 0.69 137.04 79.67 66.84 Hypo10 4.45 0.63 160.47 56.24 67.66 a root means square deviation; b squared correlation coefficient; c overall cost parameter of the PH4 pharmacophore; d cost difference between Null cost and hypothesis total cost; e lowest cost from 49 scrambled runs at a selected level of confidence of 98%, the Configuration cost = 9.85.