Exploring Novel Inhibitors for Babesia bigemina Lactate Dehydrogenase: A Computational Structural Biology Perspective

doi:10.21203/rs.3.rs-4951492/v1

Babesia bigemina is an apicomplexan parasite and causes “Texas fever” in bovines. Lactate dehydrogenase (LDH) is an essential enzyme in regulating anaerobic metabolism. The presence of five amino acid insertion in the active site of many apicomplexan LDH makes a significant difference between the parasitic LDH and the host counterpart. Therefore, apicomplexan LDH is an attractive drug target. In this study, a structure-based drug discovery approach was performed to find novel inhibitor candidates. In the first round, possible candidates were identified by following the virtual screening workflow. Then, the compounds with favourable docking scores were filtered using the QM-polarized ligand docking and induced fit docking methods. As a result, 20 novel compounds that bind to the active site of the BbigLDH and have low affinity with the host LDHs have been identified. Molecular dynamics simulations of the complexes (in total 8.8 µs) were performed and binding free energies were calculated. In conclusion, compounds named C09, C16 and C18 deserve further investigation to better understand their potential therapeutic effects on babesiosis. The findings of this study, elucidating the structural properties of BbigLDH enzyme and discovering its potential inhibitors, might pave the way for further research in developing LDH-targeted therapeutic interventions.

Lactate dehydrogenase

Babesia bigemina

drug discovery

Smith and Kilborne (1893) described the 'Texas fever' for the first time, known as babesiosis, which is caused by the apicomplexan parasite Babesia bigemina and transmitted through tick vectors [1]. The disease manifests in a spectrum of states ranging from asymptomatic carriers to severe cases marked by hemolytic anemia, fever, hemoglobinuria, and, in some instances, fatalities. Its impact extends across a diverse range of mammals, including domestic pets, farm animals, wildlife, and even humans [2], [3], [4]. B. bigemina currently prevail in tropical and subtropical regions worldwide, spanning continents such as Africa, Asia, the Americas, Southern Europe, and Australia [3]. Their presence directly influences meat and milk production, thereby influencing the competitive edge of livestock industries [5]. The primary vectors for these pathogens are Rhipicephalus (Boophilus) microplus, R. annulatus, R. australis, R. decoloratus, R. bursa, R. geigyi, and R. evertsi [2], [3]. Currently, diminazene aceturate (DA) and imidocarb dipropionate stand as the most commonly employed treatments [6], [7], [8]. However, the continued or inappropriate use of these drugs poses a significant risk of engendering drug resistance [6], [9], [10], and they also carry the potential for chemical residues in meat and milk, which can have adverse effects on human health [11], [12], [13]. This pressing concern necessitates an urgent quest for the discovery or development of alternative, effective, and economically viable drugs with minimal toxicity [14], [15].

Within apicomplexans, lactate dehydrogenase (LDH) serves as a key glycolytic enzyme, facilitating the transformation of pyruvate into lactate. [16]. It bidirectionally catalyzes the conversion of pyruvate to lactate, utilizing NADH as a co-factor and regenerating NAD + essential for glycolysis to proceed, or the conversion of lactate to pyruvate, utilizing NAD + as a co-factor and regenerating NADH [17]. The energy generated through this process is harnessed by parasites to fuel their biochemical processes and ensure their survival. The substrate specific loop, which is present in the active site, is designated for its substrate in parasitic LDH. Hosts represent canonical LDH which absence of five amino acids in the substrate specific loop differs from parasitic LDH. The five amino insertion was firstly described in LDH of Plasmodium falciparum [18]. Further studies showed that the insertion site is conserved in ancient LDH and LDH of many apicomplexan parasites including B. bigemina, but in a different amino acid sequence [5], [17], [19], [20], [21], [22], [23], [24]. The vital importance of the enzyme, and the fundamental difference between the host's LDH and the parasite's LDH in the active site makes it a more specific drug target. Further investigation of P. falciparum LDH demonstrates the potential for identifying a specific inhibitor for the enzyme [25]. Other studies related to LDH in apicomplexan parasites were published for T. gondii, C. parvum, and supports the hypothesis that apicomplexan LDH could be used as a drug target [26], [27]. The crystal structure of Babesia microti LDH (BmLDH) is the first LDH whose crystal structure has been solved in babesia genus. However, it is important to note that the BmLDH shows significant differences from apicomplexans LDHs and its structure is closely related to mammalian-like LDH [28], [29]. Yu et al., 2022 was identified a novel cDNA clone that encodes apicomplexan type LDH from Babesia orientalis (BoLDH) and solved its crystal structure in its apo form [24]. On the other hand, gossypol, an LDH inhibitor, has been explored for its impact on various Babesia species. A concentration with an IC50 value of 50 µM gossypol irreversibly halted the in vitro growth of B. bovis. Crucially, the addition of gossypol led to a noteworthy inhibition of enzymatic activity, acting as a competitive inhibitor for the binding of NADH [20]. Another study for gossypol, exhibits inhibitory impact on the in vitro growth of B. bigemina, with an IC50 of 43.97 µM after 72 hours of treatment. The maximum inhibitory concentration (IC98) is observed at 60 µM gossypol [5]. These findings demonstrate that LDH is a significant drug target for Babesias as an apicomplexan parasite.

Modern progress in computer-assisted drug development has made it possible to design drugs specifically for target proteins, thereby accelerating the creation of new medications. Structure-Based Virtual Screening has notably transformed and improved the process of drug discovery, establishing itself as one of the most effective computational techniques for drug design [30]. The principles and applications of Virtual Screening encompass the entire process, from the initial stages of receptor and library preparation, through docking and scoring, to the post-processing of the highest-scoring hits. This study focused on understanding the structure of Babesia bigemina LDH (BbigLDH) enzyme and discovering candidate inhibitors by structure-based drug discovery methods. After identifying the hit compounds using virtual screening, molecular dynamics simulations were performed to assess protein-ligand complex stability and calculate end-point binding free energies. This research could contribute to the identification of potential LDH inhibitors for the treatment of Babesia bigemina infection and to improve drug design efforts against bovine babesiosis and diseases caused by other apicomplexans.

Sequence Analysis

Protein sequences of Bos taurus LDH B (CAB96751.1), Bos taurus LDH A (BAA14170.1), Babesia bovis T2Bo (XP_001611047.1), Babesia bigemina (XP_012766822.1), Babesia ovata (XP_028866450.1), Theileria equi strain WA (XP_004833698.1), Theileria orientalis strain Shintoku (XP_009689964.1), Theileria parva strain Muguga (XP_766703.1), Theileria annulata (ADG45564.1), Plasmodium malariae (SBS88125.1), Plasmodium yoelii (XP_724101.1), Plasmodium falciparum 3D7 (XP_001349989.1), Toxoplasma gondii ME49 (XP_002368488.1) and Eimeria acervulina (ACM77785.1) were obtained from the NCBI (https://www.ncbi.nlm.nih.gov/). ClustalW (1.2.4) was used for multiple sequence alignment and the aligned sequences were visualized via ESPript 3.0 webserver [24]. The phylogenetic tree which is assessed by using the Maximum Likelihood method was constructed via MEGA 11 [31].

Homology Modelling and Model Validations

BbigLDH, whose experimentally verified 3-D structure is unknown, were generated by homology modelling using MODELLER 9.23 software [32]. Template for homology modelling was determined using NCBI BLAST tool [Basic Local Alignment Seach Tool (http://www.ncbi.nlm.nih.gov/blast)]. Resulting PDB structures with the highest max score, the lowest E-value, and the highest % identity score were evaluated. T. gondii LDH1 apo form (1PZE) [17] and ancestral apicomplexan LDH with malate (4PLC) [19] were chosen as templates for the modelling of the open and closed states of the BbigLDH enzyme, respectively. A total of 100 models were generated, and the model with the lowest DOPE (Discrete Optimized Protein Energy) and z-DOPE (normalized DOPE) scores was selected for further in silico analyses. UCSF Chimera 1.15 was employed to energy minimize the BbigLDH models, adding any missing hydrogen atoms and charges [33]. The minimization process involved 2000 steps of "Steepest Descent" followed by 1000 steps of "Conjugate Gradient," utilizing the AmberFF14SB force field.

Quality of the protein models was evaluated by using servers including ERRAT [34], ProSA [35], [36], ProCHECK / Ramanchandran Plot [37], Verify3D [38] and QMEAN [39]. To assess the reliability of specific atomic interactions, ERRAT was utilized. ProSA was employed to verify the correctness of the overall 3D homology model structure of BbigLDH and ensure proper energy optimization. Verify3D was used for 3D structure verification by computing position-specific 3D-1D scores for each residue, referencing a database of known accurate folding structures. The QMEAN tool provided a quantitative evaluation of significant geometric features in the protein structure. The Ramachandran Plot (ProCHECK) was used to visualize the φ (phi) and ψ (psi) dihedral angles of the amino acid residues, identifying conformational outliers and confirming that these angles are within the favorable regions, which are indicative of proper folding and stability.

The LDH structures of the host (Bos taurus) was obtained from the Alphafold database (Uniprod ID: (LDHA - P19858, LDHB - Q5E9B1, and LDHAL6B - Q3T056) [40] (https://alphafold.ebi.ac.uk/).

Cognate docking and Enrichment Calculations

Cognate docking was performed using the Glide XP (extra precision) scoring function and docking methodology to asses prediction quality of the docking program [41], [42]. The 3-D crystal structures of P. falciparum LDH (PfLDH) with inhibitors (PDB IDs: 1LDG [18], 1CET [43], 1T26, 1T24, 1T25 [44], 1U5A, 1U4O, 1XIV, 1U4S [45]) were obtained from the Protein Data Bank [46]. The co-crystallized ligands were split from the protein structure, then protein and ligand structures were prepared using the protein preparation workflow respectively and re-docked [47], [48]. Additionally, these ligands were analyzed for their interactions with BbigLDH using XP docking. Heavy-atom RMSD values were calculated between the docked poses and the crystal structure poses using the superposition tool from Maestro [49]. In addition, docking methodology was validated by decoys set of 1000 molecules [50]. Known 34 LDH inhibitors were obtained from the BRENDA Enzyme database [51]. 1000 decoy compounds and 34 active molecules were prepared with LigPrep [48] and docked by Glide XP mode with default settings. ROC, RIE, AUAC and BEDROC parameters were analysed by Enrichment Factor Calculation.

Structure-Based Virtual Screening, QM-Polarized Ligand Docking (QPLD) and Induced Fit Docking (IFD)

Compound library including about 570.000 compounds from the ChemDiv (Targeted diversity library, https://www.chemdiv.com/catalog/diversity-libraries/targeted-diversity-library/), Enamine (Hit Locator Library HLL-460, Discovery Diversity Set DDS-50, Discovery Diversity Set DDS-10, https://enamine.net/compound-libraries/diversity-libraries) and OTOVAchemicals (PrimScreen1, PrimScreen2, PrimScreen3, PrimScreen5, PrimScreen10 https://otavachemicals.com/products/compound-libraries-for-hts/diversity-sets) were prepared using the LigPrep [48]program. Ionization states of the compounds were predicted using the Epik at 7.5 pH [52]. Protein Preparation Workflow was performed to prepare protein structures for docking[47], [53]. Receptor grid generation was used to define the atomic positions of the binding site for each structure.

Glide-Virtual Screening Workflow was performed in three steps: High throughput virtual screening (HTVS-10%), Standard precision (SP-10%), and Extra precision (XP- the best docking scores of 100 ligands). Compounds were ranked according to XP-GScore and the best compounds having interactions with the catalytic residues and lower XP-GScores (≤-7.0 kcal/mol) were selected. Afterwards, ligands are filtered based on their ADME properties [54]. The candidates identified are within the reference values specified in the Qikprop manual. For compounds to qualify as drug candidates, several criteria and appropriate value ranges must be considered. Among these, the key five principles of Lipinski are essential [55]. Compounds with acceptable scores were then used in Quantum Mechanical Polarized Ligand Docking (QPLD). The QPLD is a computational method that combines Glide docking with QSite calculations to determine the ligand binding positions [56]. For the initial docking, the Glide XP method is employed with accurate QM level chosen for QM charges using the Jaguar module in gas phase [57], [58] Redocking is performed using the XP method. In result, ligand conformations with the lowest Emodel values obtained from the docking was selected. Compounds with XP-GScore value equal to or lower than − 7.0 kcal/mol were filtered for the induced fit docking. IFD is a molecular docking method that involves ligand repositioning onto a relaxed receptor using Prime. In this approach, a receptor grid was generated by selecting active site residues for molecular docking. Prime was utilized to perform receptor rearrangement for residues within a 5 Å distance range. The relaxed receptor was then subjected to redocking using the XP method of Glide program. The resulting poses were evaluated based on the lowest IFD scores. Compounds identified post-IFD were subjected to Glide XP docking against the structures of Bos taurus LDH A, LDH B and LDH AL6B. Compounds displaying low affinity with the host LDH structures were selected for further analysis, considering them as potential candidates for selective inhibition.

Molecular Dynamics Simulations and Binding Free Energy Calculations with MM/GBSA

MD simulation was performed for open and closed states of homology models by using the Nanoscale Molecular Dynamics (NAMD) program [59]. During the preparation of closed states, ligands were processed using the CHARMM-GUI server [60] with the inclusion of hydrogen atoms facilitated by Open Babel, followed by conversion to the mol2 format [61]. The coenzyme (NAI) was then added into the system and parametrized. The protein and ligand PDB and PSF files were generated and combined using VMD 1.9.3 software to create the complex structure [62]. The complex structure was solvated in a cubic box and subsequently neutralized, and the ionic strength was adjusted to 150 mM NaCl. Before MD simulations, 2-tier of minimization (total 5000 steps) was conducted. In the first step, constraint with constraint scaling of 100 kcal/mol/A^2 was applied to the solids, while in the second step, the entire system was minimized without any constraints. Equilibration was carried out using the NVT ensemble for 250 ps at 300 K by applying constraints (25 kcal/mol/A^2) to the solids. This was followed by an addition equilibration and relaxation step with the NPT ensemble for 275 ps at 300 K and 1 atm using the Nose–Hoover Langevin piston. The constraints on the solid were gradually reduced from 10 kcal/mol/A^2 to zero at every 25 ps. Finally, a production phase consisting of a 100 ns MD simulation was performed without any constrains using the NPT ensemble (300 K and 1 atm). All MD simulations were conducted using the CHARMM36 forcefield [63]. The MD simulation trajectories were analysed for the RMSD and RMSF values and visualized using the VMD 1.9.3 software [62].

Molecular dynamics simulations of the protein-ligand complexes were performed using the Desmond software (OPLS4 forcefield) [64] or GROMACS 2022 (CHARMM36 forcefield) [63] to asses ligand stability and to calculate binding energies of compounds. The system consisted of TIP3P water molecules in an orthorhombic box. To ensure charge neutrality and achieve an ionic strength of 150 mM NaCl, appropriate ions were added. Prior to production runs, a minimization step for 100 picoseconds and a default relaxation protocol provided by Desmond was employed to minimize and relax the system. The production runs were conducted for 20, 300, or 1000 ns at the different stages of compound selection procedure for the next stages. All production runs were performed at a temperature of 300 K and a pressure of 1 atmosphere. The integration time step was set as 2 femtoseconds. Temperature and pressure were controlled using the Nose-Hoover thermostat and Martyna-Tobias-Klein (MTK) barostat, respectively. A cutoff distance of 9 Å was applied for nonbonded interactions to calculate the short-range forces efficiently. Protein-ligand interactions were investigated with simulation interaction diagram. MM/GBSA method was employed to calculate the binding free energies (ΔG bind) using the thermal_mmgbsa.py script and the Prime program [65]. In parallel, 100 ns of MD simulations were performed using the GROMACS by applying the same parameters as above. Minimization was stopped when the maximum force < 10.0 kJ/mol. NVT equilibration (100 ps) was employed by the Particle Mesh Ewald method applied to calculate long range electrostatic interaction and the modified Berendsen thermostat (V-rescale) used to control the temperature of the simulation system. NPT equilibration (100 ps) and production run was performed using the C-rescale pressure coupling at 1 atm and V-rescale used to control the temperature of the simulation system. For the analysis of MD trajectories, GROMACS modules (gmx rms and gmx rmsf) were used, focusing on parameters such as root mean squared deviation (RMSD) and fluctuation (RMSF).

Assessment of LDH sequence similarity among apicomplexan species

Upon analysis of its nucleotide sequence, it was determined that the BbigLDH gene contains an open reading frame spanning 969 nucleotides, which encodes for a sequence of 332 amino acids. The alignment of sequences was conducted using ClustalW and visualized through ESPript 3 (Supplement file, Fig. 1). The highest sequence identity of 95.34% and 85.98% was observed between BbigLDH and Babesia ovata and Babesia bovis, respectively. Other species closely related to BbigLDH, such as Theileria orientalis exhibited a 70.66% match, while Plasmodium falciparum showed a 50.31% of identity. Notably, there is a 30% identity with Bos taurus LDH A, LDH B, which is recognized as the host for Babesia bigemina. The phylogenetic tree, constructed using the Neighbor Joining method, is provided in Supplement file, Fig. 3 [66]. The phylogenetic tree suggests that BbigLDH is positioned within the Babesia clade, with a closer evolutionary relationship to B. ovata. In contrast, LDHs originating from mammalian hosts form a distinct clade, indicating a clear separation in their evolutionary lineage from those of the apicomplexan parasites. The InterPro database analysis of BbigLDH sequence (XP_012766822.1) revealed the presence of LDH1N-NAD binding (PF0056 / 9-152) and LDH1C-alpha/beta C terminal (PF02866 / 157–315) domains. The 3D model of BbigLDH is expected to include a binding site for NADH or NAD + cofactors as well as a binding site for pyruvate or lactate substrates. The secondary structure of BbigLDH consists of sixteen alpha-helices, fourteen beta-strands, and various loops (Supplement file, Fig. 2).

BbigLDH displays a common five-amino acid insertion within its substrate specificity loop. Crystallographic analysis PfLDH have shown that this insertion contributes to a distinctive architecture in the active site, which could be a target for drug development [18]. While the sequences in this region vary across species, for BbigLDH, the sequence is DDEWT, and it is predicted to form a coil structure (Fig. 2A). BbigLDH shares identical catalytic residues with those conserved in all LDHs, which include R109, D168, H195, and R171 (residues R98, D158, H185, R161 in BbigLDH) (Supplement file, Fig. 1) [5]. Residues 93–97 (DDEWT) typically form a mobile loop in most apicomplexan LDHs (Fig. 2A). In the ternary complex, the loop moves downward to form the closed catalytic site [19]. The loop opposite to the substrate specificity loop (G236, T237, G238 in BbigLDH) is called as “opposing loop” and depicted in Fig. 2B and Supplement file, Fig. 5 [19]. The distinctive two-residue deletion within the opposing loop of modern Plasmodium LDHs, along with the shared loop structure with ancient LDH (AncLDH) and Toxoplasma LDHs, marks a significant evolutionary deviation (Supplement file, Fig. 5) [19]. Additionally, while ancestral LDH displays a slight oxaloacetate activity in common with modern Toxoplasma LDHs, this capability is absent in Plasmodium LDHs. This intriguing correlation strongly suggests that the opposing loop deletion, and potentially the presence of Ala236 and Pro246, play a pivotal role in shaping the remarkably stringent substrate specificity observed in modern Plasmodium LDHs [67]. Unlike in modern Plasmodium LDHs, a deletion in the opposing loop is not observed in BbigLDH, like Toxoplasma gondii LDH. This sets it apart from the modern Plasmodium LDHs in this aspect. In Babesia and Theileria (excluding T. equi) genus, a G[QSTN]G pattern is present (Fig. 2). Unlike in Plasmodium and B. taurus, such a pattern is not found, and there is a distinction in the opposing loop from both ancLDH and Toxoplasma LDH. While residues preceding the opposing loop in B. taurus and Plasmodium LDHs are positively charged (H and K residues), in Babesia and Theileria, these residues are observed to be aromatic (Y).

Homology modelling and model validation

The crystal structures of monomers in both closed and open states of BbigLDH were previously unknown. To address this, homology models were generated using the MODELLER 9.23 program. These models were developed for both open and closed states, incorporating pyruvate and NADH, respectively. The suitable template proteins were selected from PDB-structured proteins using BLASTp, based on the highest identity percentage. Additionally, during the selection process, emphasis was placed on ensuring that the template LDH structures belonged to apicomplexans and contained NADH and pyruvate in closed state. The template LDH proteins used for modelling included 1PZE (1.95 Å) [68] - the T. gondii LDH1 apoenzyme structure- for the open state, and 4PLC (1.50 Å) - an ancestral LDH from apicomplexans [67]- for the closed state, along with pyruvate and NADH. The generated 100 models were selected based on the lowest DOPE value and subsequently validated using various methods. One of the validation programs employed was ERRAT, which examines the consistency of atomic interactions within the modelled structure. The overall quality factors obtained from the ERRAT server for the open and closed models were 90.1274 and 93.949, respectively [69]. These values indicate a relatively high quality of the modelled structures. Another validation program utilized was ProSA [70], which evaluates the quality of generated models by constructing an energy profile of structural residue areas. The Z-scores obtained from ProSA were compared with those of experimentally validated structures to assess consistency. According to ProSA results, the Z-scores for mutant and wild-type models were − 3.3 and − 3.19, respectively. Models falling within the acceptable range of residue and energy efficiency were considered reliable. Quantitative Model Energy Analysis (QMEAN) is a valuable tool for analyzing crucial geometric aspects of protein structures, such as determining natural curvatures and calculating torsion angle potentials of amino acids [71]. The QMEAN value is originally in a range [0,1] with one being good expected from high resolution X-ray structures. In the results, QMEAN4 values were − 2.47 and − 2.19 for open state and closed state, respectively. Verify3D is a validation program that evaluates the absolute secondary and tertiary folding behaviour of a model protein by analyzing the 3D-1D scores in the context of 3D structure-sequence relationships. A higher value for amino acids is expected compared to the threshold value, and the overall simplified score determines the structure's validity [72]. Homology models passed by 88.82% and 93.79% for open state and close state, respectively according to the assessment based on favourable 3D-1D scores. A Ramachandran plot graph was generated for the model proteins using ProCHECK [73]. The percentages of residues in the favoured region were 88.1% and 86.6% respectively for the open and closed states. Based on the assessment of validation values, all these validation results confirmed that the enzymes were modelled with high accuracy and reliability (Fig. 3A). Additionally, the structural difference introduced by the 5-AA insertion in the open form of BbigLDH compared to host LDH A and B is shown in Fig. 3B. Detailed results obtained from validation were provided in Supplement file, Table 1.

Table 1

Twenty compounds, exhibiting lower affinity with BtLDH A, B and A-like B and their average values (kcal/mol)
Compound no	Compound ID	BtLDHA	BtLDHB	BtLDHA-likeB	Average
C1	Z1478445725	-4.311	-5.156	-4.875	-4.78067
C2	Z1246322905	-4.495	-5.448	-3.34	-4.42767
C3	Z826473056	-4.253	-4.534	-4.414	-4.40033
C4	Z509196536	-4.233	-5.03	-4.275	-4.51267
C5	Z1424270493	-4.559	-4.411	-5.247	-4.739
C6	Z278525320	-4.464	-4.047	-4.283	-4.26467
C7	Z1640818270	-4.863	-5.026	-4.731	-4.87333
C8	Z604033838	-3.745	-3.219	-4.344	-3.76933
C9	Z2373310039	-4.842	-5.14	-4.485	-4.82233
C10	Z973658262	-4.959	-3.858	-3.883	-4.23333
C11	Z1616793138	-4.507	-4.505	-4.672	-4.56133
C12	Z2953881032	-4.596	-3.577	-4.391	-4.188
C13	Z1253206711	-4.843	-4.798	-4.772	-4.80433
C14	Z19650139	-4.249	-4.068	-4.262	-4.193
C15	Z2400537901	-5.409	-4.512	-4.829	-4.91667
C16	Z2201943136	-4.8	-4.406	-5.435	-4.88033
C17	Z2060912694	-4.114	-4.636	-4.857	-4.53567
C18	Z2063989764	-3.427	-5.208	-5.287	-4.64067
C19	Z1531526443	-4.526	-5.325	-4.901	-4.91733
C20	Z2055833084	-4.401	-4.787	-4.808	-4.66533

Validation of Molecular Docking Procedure

The validation of the molecular docking procedure to be used in computational biology methods is an essential preliminary step to ensure the validity and reliability of the method. In the validation process, inhibitors present in the crystal structures of PfLDH with known structures were redocked into the original protein structure using the cognate docking procedure in methods section. The superposition of the resulting protein structure with the original protein structure is examined, along with the binding mode of the ligand and the RMSD values. Low RMSD values (The superimposed values for 1U4O, 1T25 and 1XIV are 0.14 Å, 0.21 Å, and 0.32 Å, respectively.) indicate a high degree of similarity between the two structures. Other RMSD values for 1U4S and 1T26 are 1.26 Å and 1.06 Å, respectively which have more higher superimposed values then other PDB structures. Additionally, when evaluating the cognate docking results, a threshold value of -7.00 kcal/mol for the docking score has been determined as the preferred threshold for assessing the HTVS docking results.

The enrichment calculation results provide valuable information about the performance of a computational model or algorithm in virtual screening or drug discovery experiments. The model was tested on a dataset consisting of 1000 ligands, including 34 known actives. The goal was to determine the model's ability to accurately rank the active compounds and enrich them at the top of the ranked list. The obtained BEDROC values ranged from 0.229 to 0.332, suggesting a moderate enrichment capability. The receiver operating characteristic (ROC) curve and its associated area under the curve (AUC) were also analyzed to evaluate the overall classification accuracy of the model. The ROC value of 0.69 suggests a reasonable discriminative ability, indicating that the model can distinguish between active and inactive compounds to a certain extent. The area under the accumulation curve (AUC) is 0.68, indicating that a significant portion of the active compounds is discovered early in the screening process. The obtained RIE value of 3.43 indicates that the model is performing better than random selection but still leaves room for improvement. It highlights the model's ability to detect structurally diverse active compounds (Supplement file, Fig. 4 and Table 2).

Structure-Based Virtual Screening, QM-Polarized Ligand Docking and Induced Fit Docking (IFD)

Four structures of BbigLDH were utilized, comprising two open and two closed states. A pair of open and closed states were generated after a 100 ns MD simulation using NAMD. For the prediction of binding sites, DoGSiteScorer method was employed on all models. [74]. In both the open and closed states of BbigLDH, druggable pockets were predicted, with a druggability score of 0.81 and simple score of 0.6. The pocket includes the pyruvate and NADH binding regions. The virtual library was initially docked to the active site of BbigLDH using GLIDE HTVS, SP and XP docking methods, respectively. The top 100 poses were then ranked based on XP GScores (Table 2). After comparison of compounds interacting with common residues, 256 ligands were selected. Subsequently, the remaining ligands underwent ADME analysis using the QikProp program. Compounds within the specified reference value range from the manual were identified (133 compounds). After ADME analysis, the Quantum Mechanics Polarized Ligand Docking method was applied, where the QM charges were calculated by using the QSite (Jaguar) and each ligand was redocked from its corresponding pose set. The poses with the lowest Glide Emodel values were chosen. Results with XP GScore values equal to or higher than − 7.00 kcal/mol were ranked to obtain the best outcomes for 70 models. Following this, induced fit docking was employed This involved relaxing the ligand into a conformationally relaxed receptor state, followed by redocking using Glide. The resulting 70 poses were ranked based on their IFD score, with lower scores being prioritized. Remaining compounds were also docked to the host LDH (Bos taurus LDH A, LDH B and LDH AL6B) to identify selective inhibitor candidates. The structural difference introduced by the 5-AA insertion in the open form of BbigLDH compared to host LDH A and B is shown in Fig. 3B. Twenty compounds, exhibiting lower affinity with host LDHs and available for purchase were selected for molecular dynamics simulations and end-point binding free energy calculations (Table 1).

The evaluation of potential drug candidates involves a thorough assessment of various criteria, with a close examination of their physicochemical properties intricately linked to their structural attributes. Minor modifications in their structures have resulted in enhanced properties. In our study, the hit compounds underwent analysis to determine their compliance with Lipinski's rule of five. [55]. According to Lipinski's rule of five, an established set of guidelines for drug likeness, the molecular weight of potential drug candidates should be less than 500 Da. Additionally, the hydrogen bond donor count should be fewer than 5, the hydrogen bond acceptor count should be less than 10, and the predicted octanol/water partition coefficient should ideally fall within the range of -2.0 to 4. These criteria help assess the likelihood of a compound's success as a drug candidate based on its physicochemical properties [75]. All values for the compounds fall within the specified range for Lipinski rules and PlogP_0/w. The Plog values for both hERG channel blocking and serum albumin binding fell within the reference ranges [76]. The predicted IC50 value for blocking HERG K + channels should fall within the accepted range, with values greater than − 5.0. In this study, the compound with the largest value is C18 (5.073), while the one with the smallest value is compound C9 (-2.504). On the other hand, the predicted apparent Caco-2 cell permeability should be assessed, with values indicating low permeability if less than 25 and high permeability if greater than 500 [77]. According to this information, it is observed that the compound with the highest permeability is C20 (2104.635), while the one with the least permeable is C19 (74.968). The predicted blood-brain barrier permeability coefficient should be within the range of -3.0 to 1.2 [78]. The predicted binding capacity to human serum albumin is considered appropriate if it falls between − 1.5 and 1.5 [79]. All values for the compounds fall within the specified range for PlogBB, PlogKhsa. The predicted human oral absorption, rated on a scale of 0 to 100, is considered favorable if it exceeds 70%. The compound with the highest human oral absorption rates is C18 and C20 (both 100), whereas C19 has a rate of 72.307 (Supplement file, Table 3).

Table 2

List of compounds which have favourable XP GScores. Free binding energy (∆G) values were listed for the last 5 ns of 20 ns MD simulation.
Compound no	Compound ID	Structure of BbigLDH	HTVS; XP gscore	QPLD; XP gscore	IFD; XP gscore	Prime energy	IFD score	dG Average
C1	Z1478445725	Closed state after MDS	-8.343	-8.515	-9.569	-13590.36	-689.09	-67.3425
C2	Z1246322905	Closed state after MDS	-8.608	-8.140	-9.048	-13502.23	-684.16	-58.4571
C3	Z826473056	Closed state after MDS	-8.236	-8.530	-11.017	-13460.85	-684.06	-64.0088
C4	Z509196536	Closed state	-8.997	-8.389	-10.565	-13452.59	-683.19	-57.2994
C5	Z1424270493	Closed state	-8.869	-7.598	-9.404	-13471.78	-682.99	-60.7591
C6	Z278525320	Closed state after MDS	-9.497	-7.480	-10.557	-13439.77	-682.55	-54.8145
C7	Z1640818270	Closed state after MDS	-8.479	-8.473	-7.652	-13462.95	-680.8	-64.4233
C8	Z604033838	Closed state	-8.382	-8.250	-10.285	-13409.02	-680.74	-66.175
C9	Z2373310039	Closed state after MDS	-8.173	-7.958	-7.998	-13442.56	-680.13	-75.825
C10	Z973658262	Closed state	-9.176	-9.182	-9.803	-13401.97	-679.9	-46.3542
C11	Z1616793138	Closed state	-8.648	-8.555	-9.624	-13403.46	-679.8	-46.0935
C12	Z2953881032	Closed state after MDS	-8.202	-7.972	-8.058	-13428.81	-679.5	-71.5432
C13	Z1253206711	Closed state	-8.925	-7.731	-10.617	-13374.26	-679.33	-53.7573
C14	Z19650139	Closed state	-8.701	-9.111	-7.606	-13429.61	-679.09	-60.4939
C15	Z2400537901	Closed state	-10.464	-8.919	-9.831	-13379.26	-678.79	-67.0062
C16	Z2201943136	Open state	-8.992	-8.443	-8.72	-13399.65	-678.70	-60.9869
C17	Z2060912694	Open state	-8.774	-8.274	-10.44	-13364.38	-678.66	-60.6208
C18	Z2063989764	Open state	-9.105	-7.750	-9.56	-13377.59	-678.43	-72.3022
C19	Z1531526443	Open state	-8.641	-8.706	-9.55	-13338.36	-676.47	-72.2696
C20	Z2055833084	Open state	-8.482	-7.328	-9.19	-13311.35	-674.76	-44.1338

Compound prioritization based on molecular dynamics simulations and binding energy

In this study, we conducted MD simulations and binding free energy calculations of 20 protein-ligand complexes obtained from the docking campaigns using the Desmond program for 20 ns. Furthermore, to determine ligand stability over a long period of time, MDS with GROMACS was performed for 100 ns, too.

Based on the RMSD values from the GROMACS simulations and the ∆G average values from the 20 ns MD simulations using Desmond (Table 2), we identified C7, C9, C12, C14, C15, C16, C18 and C20 compounds as candidates for an extended MD simulation for 300 ns. Subsequently, we assessed the RMSD values obtained from a 300 ns MD simulation (Fig. 4). According to the Lig-fit-prot RMSD values in Fig. 4B, all compounds except C20 were found in the binding cavity at the end of the simulation. Moreover, C16 demonstrated remarkable positional stability, maintaining a 3 Å RMSD value throughout the simulation. According to the MM/GBSA ∆G Distribution graph in Fig. 6A, both C9 and subsequently C16 displayed favourable binding energies at the last 30 ns of the 300 ns simulation.

In the subsequent phase, we selected C9, C15, C16, and C18 for an extensive 1000 ns MDS analysis according to the ∆G and RMSD values, and interaction with the active site residues. The RMSD analysis graph in Fig. 5B reveals that C16 maintained stability within the range of 2–4 Å. In the case of C9, stability was observed after 400 ns, with a stability range of 4–6 Å. C18 exhibited stability after 150 ns, within a range of 6–8 Å. C15 moving away from the binding pocket, therefore RMSD analysis was not performed. In Fig. 5C, the RMSF (Root Mean Square Fluctuation) graph indicates that residues between 80 and 100, which include the substrate-specific loop, exhibit dramatic changes for C18 compared to other compounds. Additionally, the residues between 235 and 240, which encompass the opposing loop, show changes for C9 and C16, distinct from those observed for C18.

According to the MM/GBSA ∆G Distribution graph in Fig. 6B, both C9, C16 and C18 displayed favourable binding energies in the last 100 ns of the 1000 ns simulation. The ∆G averages are − 52.5 kcal/mol, -62.1 kcal/mol and − 62.5 kcal/mol for C9, C16 and C18, respectively.

Figure 7 illustrates the average energies contributing to the overall MM-GBSA binding energy. Calculations were based on snapshots taken from the last 100 ns of a 1 µs MD simulation. The MM-GBSA ΔG binding energies are quite similar for C16, and C18 (− 62.11 kcal/mol, and − 62.52 kcal/mol, respectively), but the main driving forces behind these binding energies different. For C9-protein binding, lipophilic energy and van der Waals (vdW) energy are the predominant types of energy. In contrast, the contributions of coulomb and hydrogen bond (H-bond) energies in the binding of C16 and C18 are markedly higher from those in C9-protein binding.

In Fig. 8, all three compounds are found in the expected regions after the 1000 ns simulation time. All compounds are positioned within the enzyme's catalytic pocket, which is of crucial importance, and interactions with catalytic residues are evident. C9 retained its significant H-bond interactions after IFD, but lost them after MD. In Supplementary Fig. 6A and 7A, according to BbigLDH-C9 interactions, there are water bridges which interaction mediated between a water molecule and residues including T82, L85, D93, I128, and T129. On the other hand, upon examination of Supplement file, Fig. 6B and 7B, it is evident that in the BbigLDH-C16 complex, there are hydrogen bond interactions with residues N17 (77%), L85 (97%), R87 (94%), and N130 (70%), along with water Bridges involving L231 and pi-pi stacking with Y241 (32%). In Supplement file, Fig. 6C and 7C, hydrogen bond interactions are observed with residues R87, I128, N130 (46%) and R161 (39%) in the BbigLDH-C18 complex, along with pi-cation interactions involving R161 (33%). In Fig. 9A, after the induced fit docking stage, interactions with R87, I128, and N130 residues were generally observed with all compounds. On the other hand in Fig. 9B, after the 1000 ns simulation times, interactions with R87, N130 conserved. It can be observed that the interaction with these residues is maintained throughout the MD simulation for both C16 and C18 (Supplement file, Fig. 6). It is believed that these residues play a significant role in inhibition against BbigLDH. In the cognate docking stage, inhibitors found in the open and closed PfLDH structures, whose PDB structures are known, were docked to the open and closed states of BbigLDH structures. The interaction of the original PDB structures with the inhibitor and the interaction of the BbigLDH structures with the same inhibitor are shown in Supplement file, Fig. 8. It is observed that the lack of a stable region in the substrate-specific loop region of the open state of PfLDH structures results in differences with the open state of BbigLDH structure. Additionally, it is seen that the residues interacting with the ligands in the open state of BbigLDH structure and also known to have catalytic properties in some cases are R87, I128, N130, R161, and S239. On the other hand, for the closed state of BbigLDH, the most important residues interacting with the ligands are R98, N130, R161, H185, S239, and H240, as shown in Supplement file, Fig. 9. In conclusion, it is considered that compounds showing similarity in the critical residues obtained from the redocking of PfLDH inhibitors, whose PDB structure is known, to BbigLDH, and interacting with residues C09, C16, and C18, may have similar inhibitory properties on BbigLDH.

These findings strongly indicate that compounds located in the pocket where the cofactor and substrate bind may inhibit enzyme activity through these specific residue interactions. Particularly, C16 may be a promising candidate for in vitro inhibitor studies, given its stable interaction with BbigLDH. Compounds Z2373310039 (C09), Z2201943136 (C16) and Z2063989764 (C18) may warrant further investigation to comprehend their potential therapeutic effects on babesiosis.

Recent reports have indicated the emergence of resistance against certain existing anti-babesial drugs. This highlights the urgency of identifying targets for the development of new drugs with distinct modes of action. The lactate dehydrogenase enzyme plays a pivotal role in the anaerobic pathway of the parasite, making it a promising candidate for further drug design endeavours. We employed homology modelling to predict the unknown structure of BbigLDH, followed by structure-based virtual screening. Docking scores, interaction with active site residues, affinity to host LDHs, ADME properties, and stability through simulations were considered during identifying of potential inhibitors. The stability of the three candidate compounds (C9, C16, and C18) was validated through long molecular dynamics simulations and further binding free energy calculations. The obtained results provided insights into the molecular foundation of interactions and, consequently, the mechanism of inhibition. Taken together, the comprehensive analysis of these findings underscores the utility of structure-based virtual screening in providing crucial insights into the essential structural and binding features necessary for the design of novel BbigLDH inhibitors, offering potential avenues for the treatment of bovine babesiosis.

Ethical Approval

Not applicable.

Funding

This work has been supported by Marmara University-Scientific Research Projects Coordination Unit under grant number: FDK-2022-10061.

Author contributions

Safiye Merve Bostancioglu: Investigation; Visualization; Writing-Original Draft.

Ozal Mutlu: Conceptualization; Methodology; Investigation; Visualization; Writing-Original Draft.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Disclosure statement

No potential conflict of interest was reported by the authors.

Availability of data and materials

No datasets were generated or analysed during the current study.

ORCID

Safiye Merve Bostancioglu

https://orcid.org/0000-0002-0239-3966

Ozal Mutlu

https://orcid.org/0000-0003-4551-5780

T. Smith and F. L. Kilborne, Investigations into the nature, causation, and prevention of Texas or southern cattle fever /. Washington, D.C. : U.S. Dept. of Agriculture, Bureau of Animal Industry, 1893. doi: 10.5962/bhl.title.124068.
L. Schnittger, A. E. Rodriguez, M. Florin-Christensen, and D. A. Morrison, ‘Babesia: A world emerging’, Infect. Genet. Evol., vol. 12, no. 8, pp. 1788–1809, Dec. 2012, doi: 10.1016/j.meegid.2012.07.004.
R. Bock, L. Jackson, A. De Vos, and W. Jorgensen, ‘Babesiosis of cattle’, Parasitology, vol. 129, no. S1, pp. S247–S269, Oct. 2004, doi: 10.1017/S0031182004005190.
M. Calvopiña et al., ‘Case report: First report on human infection by tick-borne Babesia bigemina in the Amazon region of Ecuador’, Front. Public Health, vol. 11, p. 1079042, Aug. 2023, doi: 10.3389/fpubh.2023.1079042.
L. He, R. G. Bastos, L. Yu, J. M. Laughery, and C. E. Suarez, ‘Lactate Dehydrogenase as a Potential Therapeutic Drug Target to Control Babesia bigemina’, Front. Cell. Infect. Microbiol., vol. 12, p. 870852, Apr. 2022, doi: 10.3389/fcimb.2022.870852.
J. Mosqueda, A. Olvera-Ramirez, G. Aguilar-Tipacamu, and G. J. Canto, ‘Current Advances in Detection and Treatment of Babesiosis’, Curr. Med. Chem., vol. 19, no. 10, pp. 1504–1518, Apr. 2012, doi: 10.2174/092986712799828355.
J. H. D. Silva et al., ‘Chemoprophylaxis for babesiosis and anaplasmosis in cattle: case report’, Rev. Bras. Parasitol. Veterinária, vol. 29, no. 4, p. e010520, 2020, doi: 10.1590/s1984-29612020096.
N. D. Shah, A. U. Bhikane, R. K. Jadhav, S. G. Chavhan, and A. Mohan, ‘Therapeutic management of Babesiosis alone and its mixed infection with Theileriosis in Cattle’, 2019.
B. Tuvshintulga, T. Sivakumar, N. Yokoyama, and I. Igarashi, ‘Development of unstable resistance to diminazene aceturate in Babesia bovis’, Int. J. Parasitol. Drugs Drug Resist., vol. 9, pp. 87–92, Apr. 2019, doi: 10.1016/j.ijpddr.2019.02.001.
Todorovic R. A., ‘Chemourophylaxis (Imidocarb. Against Babesia bigemina and Babesia argentina Infections’, Am J Vet Res, vol. 34(9):1153-61, Sep. 1973, doi: PMID: 4747036.
N. G. Coldham et al., ‘Imidocarb residues in edible bovine tissues and in vitro assessment of imidocarb metabolism and cytotoxicity’, Drug Metab. Dispos. Biol. Fate Chem., vol. 23, no. 4, pp. 501–505, Apr. 1995.
H. J. Vial and A. Gorenflot, ‘Chemotherapy against babesiosis’, Vet. Parasitol., vol. 138, no. 1–2, pp. 147–160, May 2006, doi: 10.1016/j.vetpar.2006.01.048.
A. Zintl, G. Mulcahy, H. E. Skerrett, S. M. Taylor, and J. S. Gray, ‘Babesia divergens , a Bovine Blood Parasite of Veterinary and Zoonotic Importance’, Clin. Microbiol. Rev., vol. 16, no. 4, pp. 622–636, Oct. 2003, doi: 10.1128/CMR.16.4.622-636.2003.
C. E. Suarez and S. Noh, ‘Emerging perspectives in the research of bovine babesiosis and anaplasmosis’, Vet. Parasitol., vol. 180, no. 1–2, pp. 109–125, Aug. 2011, doi: 10.1016/j.vetpar.2011.05.032.
N. McHARDY, R. M. Woollon, R. B. Clampitt, J. A. James, and R. J. Crawley, ‘Efficacy, toxicity and metabolism of imidocarb dipropionate in the treatment of Babesia ovis infection in sheep’, Res. Vet. Sci., vol. 41, no. 1, pp. 14–20, Jul. 1986, doi: 10.1016/S0034-5288(18)30565-4.
C. A. Homewood and K. D. Neame, ‘Conversion of glucose to lactate by intraerythrocytic Plasmodium berghei’, Ann. Trop. Med. Parasitol., vol. 77, no. 2, pp. 127–129, Apr. 1983, doi: 10.1080/00034983.1983.11811686.
K. L. Kavanagh, R. A. Elling, and D. K. Wilson, ‘Structure of Toxoplasma gondii LDH1: Active-Site Differences from Human Lactate Dehydrogenases and the Structural Basis for Efficient APAD ⁺ Use ^,’ Biochemistry, vol. 43, no. 4, pp. 879–889, Feb. 2004, doi: 10.1021/bi035108g.
C. R. Dunn et al., ‘The structure of lactate dehydrogenase from Plasmodium falciparum reveals a new target for anti-malarial design’, Nat. Struct. Biol., vol. 3, no. 11, pp. 912–915, Nov. 1996, doi: 10.1038/nsb1196-912.
J. I. Boucher, J. R. Jacobowitz, B. C. Beckett, S. Classen, and D. L. Theobald, ‘An atomic-resolution view of neofunctionalization in the evolution of apicomplexan lactate dehydrogenases’, eLife, vol. 3, p. e02304, Jun. 2014, doi: 10.7554/eLife.02304.
S. Bork, M. Okamura, S. Boonchit, H. Hirata, N. Yokoyama, and I. Igarashi, ‘Identification of Babesia bovis l-lactate dehydrogenase as a potential chemotherapeutical target against bovine babesiosis’, Mol. Biochem. Parasitol., vol. 136, no. 2, pp. 165–172, Aug. 2004, doi: 10.1016/j.molbiopara.2004.03.009.
D. Turgut-Balik et al., ‘Cloning, sequence and expression of the lactate dehydrogenase gene from the human malaria parasite, Plasmodium vivax’, Biotechnol. Lett., vol. 26, no. 13, pp. 1051–1055, Jul. 2004, doi: 10.1023/B:BILE.0000032958.78158.10.
O. Mutlu and D. Turgut Balik, ‘Kinetic Analysis of the Amino Terminal End of Active Site Loop of Lactate Deyhdrogenase from Plasmodium Vivax’, Balk. Med. J., vol. 29, no. 4, pp. 364–369, Dec. 2012, doi: 10.5152/balkanmedj.2012.047.
A. Erdemir and O. Mutlu, ‘Functional and structural characterization of the pentapeptide insertion of Theileria annulata lactate dehydrogenase by site-directed mutagenesis, comparative modeling and molecular dynamics simulations’, J. Mol. Graph. Model., vol. 74, pp. 344–351, Jun. 2017, doi: 10.1016/j.jmgm.2017.04.010.
L. Yu, Q. Liu, W. Luo, J. Zhao, H. F. Alzan, and L. He, ‘The Structural Basis of Babesia orientalis Lactate Dehydrogenase’, Front. Cell. Infect. Microbiol., vol. 11, p. 790101, Jan. 2022, doi: 10.3389/fcimb.2021.790101.
M. S. Gomez et al., ‘Substrate and cofactor specificity and selective inhibition of lactate dehydrogenase from the malarial parasite P. falciparum’, Mol. Biochem. Parasitol., vol. 90, no. 1, pp. 235–246, Dec. 1997, doi: 10.1016/S0166-6851(97)00140-0.
F. Al-Anouti, S. Tomavo, S. Parmley, and S. Ananvoranich, ‘The Expression of Lactate Dehydrogenase Is Important for the Cell Cycle of Toxoplasma gondii’, J. Biol. Chem., vol. 279, no. 50, pp. 52300–52311, Dec. 2004, doi: 10.1074/jbc.M409175200.
D. Madern, ‘Evolution of Cryptosporidium parvum Lactate Dehydrogenase from Malate Dehydrogenase by a Very Recent Event of Gene Duplication’, Mol. Biol. Evol., vol. 21, no. 3, pp. 489–497, Dec. 2003, doi: 10.1093/molbev/msh042.
L. Yu et al., ‘Crystal structures of Babesia microti lactate dehydrogenase BmLDH reveal a critical role for Arg99 in catalysis’, FASEB J., vol. 33, no. 12, pp. 13669–13682, Dec. 2019, doi: 10.1096/fj.201901259R.
S. Akash et al., ‘Anti-parasitic drug discovery against Babesia microti by natural compounds: an extensive computational drug design approach’, Front. Cell. Infect. Microbiol., vol. 13, p. 1222913, Aug. 2023, doi: 10.3389/fcimb.2023.1222913.
E. Lionta, G. Spyrou, D. Vassilatis, and Z. Cournia, ‘Structure-Based Virtual Screening for Drug Discovery: Principles, Applications and Recent Advances’, Curr. Top. Med. Chem., vol. 14, no. 16, pp. 1923–1938, Oct. 2014, doi: 10.2174/1568026614666140929124445.
K. Tamura, G. Stecher, and S. Kumar, ‘MEGA11: Molecular Evolutionary Genetics Analysis Version 11’, Mol. Biol. Evol., vol. 38, no. 7, pp. 3022–3027, Jun. 2021, doi: 10.1093/molbev/msab120.
B. Webb and A. Sali, ‘Comparative Protein Structure Modeling Using MODELLER’, Curr. Protoc. Bioinforma., vol. 54, no. 1, Jun. 2016, doi: 10.1002/cpbi.3.
E. F. Pettersen et al., ‘UCSF Chimera—A visualization system for exploratory research and analysis’, J. Comput. Chem., vol. 25, no. 13, pp. 1605–1612, Oct. 2004, doi: 10.1002/jcc.20084.
C. Colovos and T. O. Yeates, ‘Verification of protein structures: Patterns of nonbonded atomic interactions’, Protein Sci., vol. 2, no. 9, pp. 1511–1519, Sep. 1993, doi: 10.1002/pro.5560020916.
M. J. Sippl, ‘Recognition of errors in three-dimensional structures of proteins’, Proteins Struct. Funct. Genet., vol. 17, no. 4, pp. 355–362, Dec. 1993, doi: 10.1002/prot.340170404.
M. Wiederstein and M. J. Sippl, ‘ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins’, Nucleic Acids Res., vol. 35, no. Web Server, pp. W407–W410, May 2007, doi: 10.1093/nar/gkm290.
R. A. Laskowski, M. W. MacArthur, D. S. Moss, and J. M. Thornton, ‘PROCHECK: a program to check the stereochemical quality of protein structures’, J. Appl. Crystallogr., vol. 26, no. 2, pp. 283–291, Apr. 1993, doi: 10.1107/S0021889892009944.
J. U. Bowie, R. Lüthy, and D. Eisenberg, ‘A Method to Identify Protein Sequences That Fold into a Known Three-Dimensional Structure’, Science, vol. 253, no. 5016, pp. 164–170, Jul. 1991, doi: 10.1126/science.1853201.
P. Benkert, S. C. E. Tosatto, and D. Schomburg, ‘QMEAN: A comprehensive scoring function for model quality assessment’, Proteins Struct. Funct. Bioinforma., vol. 71, no. 1, pp. 261–277, Apr. 2008, doi: 10.1002/prot.21715.
J. Jumper et al., ‘Highly accurate protein structure prediction with AlphaFold’, Nature, vol. 596, no. 7873, pp. 583–589, Aug. 2021, doi: 10.1038/s41586-021-03819-2.
R. A. Friesner et al., ‘Extra Precision Glide: Docking and Scoring Incorporating a Model of Hydrophobic Enclosure for Protein−Ligand Complexes’, J. Med. Chem., vol. 49, no. 21, pp. 6177–6196, Oct. 2006, doi: 10.1021/jm051256o.
Schrödinger Release 2022-4: Glide, Schrödinger, LLC, New York, NY, 2022.
J. A. Read, K. W. Wilkinson, R. Tranter, R. B. Sessions, and R. L. Brady, ‘Chloroquine Binds in the Cofactor Binding Site ofPlasmodium falciparum Lactate Dehydrogenase’, J. Biol. Chem., vol. 274, no. 15, pp. 10213–10218, Apr. 1999, doi: 10.1074/jbc.274.15.10213.
A. Cameron et al., ‘Identification and Activity of a Series of Azole-based Compounds with Lactate Dehydrogenase-directed Anti-malarial Activity’, J. Biol. Chem., vol. 279, no. 30, pp. 31429–31439, Jul. 2004, doi: 10.1074/jbc.M402433200.
R. Conners et al., ‘Mapping the binding site for gossypol-like inhibitors of Plasmodium falciparum lactate dehydrogenase’, Mol. Biochem. Parasitol., vol. 142, no. 2, pp. 137–148, Aug. 2005, doi: 10.1016/j.molbiopara.2005.03.015.
H. M. Berman, ‘The Protein Data Bank’, Nucleic Acids Res., vol. 28, no. 1, pp. 235–242, Jan. 2000, doi: 10.1093/nar/28.1.235.
Schrödinger Release 2022-4: Protein Preparation Wizard; Epik, Schrödinger, LLC, New York, NY, 2022; Impact, Schrödinger, LLC, New York, NY; Prime, Schrödinger, LLC, New York, NY, 2022.
Schrödinger Release 2022-4: LigPrep, Schrödinger, LLC, New York, NY, 2022.
Schrödinger Release 2022-4: Maestro, Schrödinger, LLC, New York, NY, 2022.
R. A. Friesner et al., ‘Glide: A New Approach for Rapid, Accurate Docking and Scoring. 1. Method and Assessment of Docking Accuracy’, J. Med. Chem., vol. 47, no. 7, pp. 1739–1749, Mar. 2004, doi: 10.1021/jm0306430.
A. Chang et al., ‘BRENDA, the ELIXIR core data resource in 2021: new developments and updates’, Nucleic Acids Res., vol. 49, no. D1, pp. D498–D508, Jan. 2021, doi: 10.1093/nar/gkaa1025.
J. C. Shelley, A. Cholleti, L. L. Frye, J. R. Greenwood, M. R. Timlin, and M. Uchimaya, ‘Epik: a software program for pK a prediction and protonation state generation for drug-like molecules’, J. Comput. Aided Mol. Des., vol. 21, no. 12, pp. 681–691, Dec. 2007, doi: 10.1007/s10822-007-9133-z.
G. Madhavi Sastry, M. Adzhigirey, T. Day, R. Annabhimoju, and W. Sherman, ‘Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments’, J. Comput. Aided Mol. Des., vol. 27, no. 3, pp. 221–234, Mar. 2013, doi: 10.1007/s10822-013-9644-8.
Schrödinger Release 2022-4: QikProp, Schrödinger, LLC, New York, NY, 2022-4.
C. A. Lipinski, ‘Drug-like properties and the causes of poor solubility and poor permeability’, J. Pharmacol. Toxicol. Methods, vol. 44, no. 1, pp. 235–249, Jul. 2000, doi: 10.1016/S1056-8719(00)00107-6.
A. E. Cho, V. Guallar, B. J. Berne, and R. Friesner, ‘Importance of accurate charges in molecular docking: Quantum mechanical/molecular mechanical (QM/MM) approach’, J. Comput. Chem., vol. 26, no. 9, pp. 915–931, Jul. 2005, doi: 10.1002/jcc.20222.
A. D. Bochevarov et al., ‘Jaguar: A high‐performance quantum chemistry software program with strengths in life and materials sciences’, Int. J. Quantum Chem., vol. 113, no. 18, pp. 2110–2142, Sep. 2013, doi: 10.1002/qua.24481.
Schrödinger Release 2022-4: Jaguar, Schrödinger, LLC, New York, NY, 2022.
J. C. Phillips et al., ‘Scalable molecular dynamics on CPU and GPU architectures with NAMD’, J. Chem. Phys., vol. 153, no. 4, p. 044130, Jul. 2020, doi: 10.1063/5.0014475.
J. Lee et al., ‘CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field’, J. Chem. Theory Comput., vol. 12, no. 1, pp. 405–413, Jan. 2016, doi: 10.1021/acs.jctc.5b00935.
N. M. O’Boyle, M. Banck, C. A. James, C. Morley, T. Vandermeersch, and G. R. Hutchison, ‘Open Babel: An open chemical toolbox’, J. Cheminformatics, vol. 3, no. 1, p. 33, Dec. 2011, doi: 10.1186/1758-2946-3-33.
W. Humphrey, A. Dalke, and K. Schulten, ‘VMD: Visual molecular dynamics’, J. Mol. Graph., vol. 14, no. 1, pp. 33–38, Feb. 1996, doi: 10.1016/0263-7855(96)00018-5.
K. Vanommeslaeghe et al., ‘CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields’, J. Comput. Chem., p. NA-NA, 2009, doi: 10.1002/jcc.21367.
Schrödinger Release 2022-4: Desmond Molecular Dynamics System, D. E. Shaw Research, New York, NY, 2022.
M. P. Jacobson, R. A. Friesner, Z. Xiang, and B. Honig, ‘On the Role of the Crystal Environment in Determining Protein Side-chain Conformations’, J. Mol. Biol., vol. 320, no. 3, pp. 597–608, Jul. 2002, doi: 10.1016/S0022-2836(02)00470-9.
Naruya Saitou and M. Nei, ‘The neighbor-joining method: a new method for reconstructing phylogenetic trees.’, Mol. Biol. Evol., Jul. 1987, doi: 10.1093/oxfordjournals.molbev.a040454.
J. I. Boucher, J. R. Jacobowitz, B. C. Beckett, S. Classen, and D. L. Theobald, ‘An atomic-resolution view of neofunctionalization in the evolution of apicomplexan lactate dehydrogenases’, eLife, vol. 3, p. e02304, Jun. 2014, doi: 10.7554/eLife.02304.
K. L. Kavanagh, R. A. Elling, and D. K. Wilson, ‘Structure of Toxoplasma gondii LDH1: Active-Site Differences from Human Lactate Dehydrogenases and the Structural Basis for Efficient APAD ⁺ Use ^,’ Biochemistry, vol. 43, no. 4, pp. 879–889, Feb. 2004, doi: 10.1021/bi035108g.
C. Colovos and T. O. Yeates, ‘Verification of protein structures: patterns of nonbonded atomic interactions’, Protein Sci. Publ. Protein Soc., vol. 2, no. 9, pp. 1511–1519, Sep. 1993, doi: 10.1002/pro.5560020916.
M. Wiederstein and M. J. Sippl, ‘ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins’, Nucleic Acids Res., vol. 35, no. Web Server issue, pp. W407-410, Jul. 2007, doi: 10.1093/nar/gkm290.
P. Benkert, S. C. E. Tosatto, and D. Schomburg, ‘QMEAN: A comprehensive scoring function for model quality assessment’, Proteins, vol. 71, no. 1, pp. 261–277, Apr. 2008, doi: 10.1002/prot.21715.
J. U. Bowie, R. Lüthy, and D. Eisenberg, ‘A method to identify protein sequences that fold into a known three-dimensional structure’, Science, vol. 253, no. 5016, pp. 164–170, Jul. 1991, doi: 10.1126/science.1853201.
R. A. Laskowski, M. W. MacArthur, D. S. Moss, and J. M. Thornton, ‘PROCHECK: a program to check the stereochemical quality of protein structures’, J. Appl. Crystallogr., vol. 26, no. 2, pp. 283–291, Apr. 1993, doi: 10.1107/S0021889892009944.
A. Volkamer, D. Kuhn, F. Rippmann, and M. Rarey, ‘DoGSiteScorer: a web server for automatic binding site prediction, analysis and druggability assessment’, Bioinformatics, vol. 28, no. 15, pp. 2074–2075, Aug. 2012, doi: 10.1093/bioinformatics/bts310.
W. L. Jorgensen and E. M. Duffy, ‘Prediction of drug solubility from Monte Carlo simulations’, Bioorg. Med. Chem. Lett., vol. 10, no. 11, pp. 1155–1158, Jun. 2000, doi: 10.1016/S0960-894X(00)00172-4.
A. Cavalli, E. Poluzzi, F. De Ponti, and M. Recanatini, ‘Toward a Pharmacophore for Drugs Inducing the Long QT Syndrome: Insights from a CoMFA Study of HERG K ⁺ Channel Blockers’, J. Med. Chem., vol. 45, no. 18, pp. 3844–3853, Aug. 2002, doi: 10.1021/jm0208875.
M. Yazdanian, S. L. Glynn, J. L. Wright, and A. Hawi, ‘Correlating Partitioning and Caco-2 Cell Permeability of Structurally Diverse Small Molecular Weight Compounds’, Pharm. Res., vol. 15, no. 9, pp. 1490–1494, 1998, doi: 10.1023/A:1011930411574.
J. M. Luco, ‘Prediction of the Brain−Blood Distribution of a Large Set of Drugs from Structurally Derived Descriptors Using Partial Least-Squares (PLS) Modeling’, J. Chem. Inf. Comput. Sci., vol. 39, no. 2, pp. 396–404, Mar. 1999, doi: 10.1021/ci980411n.
G. Colmenarejo, A. Alvarez-Pedraglio, and J.-L. Lavandera, ‘Cheminformatic Models To Predict Binding Affinities to Human Serum Albumin’, J. Med. Chem., vol. 44, no. 25, pp. 4370–4378, Dec. 2001, doi: 10.1021/jm010960b.

No competing interests reported.

Supplementaryfiles.docx

Exploring Novel Inhibitors for Babesia bigemina Lactate Dehydrogenase: A Computational Structural Biology Perspective

Status:

Version 1

Abstract

Figures

Introduction

Materials & Methods

Sequence Analysis

Homology Modelling and Model Validations

Cognate docking and Enrichment Calculations

Structure-Based Virtual Screening, QM-Polarized Ligand Docking (QPLD) and Induced Fit Docking (IFD)

Molecular Dynamics Simulations and Binding Free Energy Calculations with MM/GBSA

Results & Discussions

Assessment of LDH sequence similarity among apicomplexan species

Homology modelling and model validation

Validation of Molecular Docking Procedure

Structure-Based Virtual Screening, QM-Polarized Ligand Docking and Induced Fit Docking (IFD)

Compound prioritization based on molecular dynamics simulations and binding energy

Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1