key: cord-0065788-sql3ce3t
authors: Hema, Kanipakam; Ahamad, Shahzaib; Joon, Hemant Kumar; Pandey, Rajan; Gupta, Dinesh
title: Atomic Resolution Homology Models and Molecular Dynamics Simulations of Plasmodium falciparum Tubulins
date: 2021-06-30
journal: ACS Omega
DOI: 10.1021/acsomega.1c01988
sha: 5c3644b5ce8ad196a7911f9f0178ac84fd3d695e
doc_id: 65788
cord_uid: sql3ce3t

[Image: see text] Microtubules are tubulin polymers present in the eukaryotic cytoskeleton essential for structural stability and cell division that are also roadways for intracellular transport of vesicles and organelles. In the human malaria parasite Plasmodium falciparum, apart from providing structural stability and cell division, microtubules also facilitate important biological activities crucial for parasite survival in hosts, such as egression and motility. Hence, parasite structures and processes involving microtubules are among the most important drug targets for discovering much-needed novel Plasmodium inhibitors. The current study aims to construct reliable and high-quality 3D models of α-, β-, and γ-tubulins using various modeling techniques. We identified a common binding pocket specific to Plasmodium α-, β-, and γ-tubulins. Molecular dynamics simulations confirmed the stability of the Plasmodium tubulin 3D structures. The models generated in the present study may be used for protein–protein and protein–drug interaction investigations targeted toward designing malaria parasite tubulin-specific inhibitors.

Malaria is one of the most devastating global health problems due to the increased morbidity and mortality rates. Globally, there were an estimated 229 million malaria cases and 409,000 deaths reported by WHO (2020). 1 Among malaria-causing parasites, Plasmodium falciparum (P. falciparum) causes the most virulent form of human malaria. The risk of severe illnesses is higher in children and pregnant women in the African zone, South-East Asia, and Eastern Mediterranean zones. 1, 2 Plasmodium parasites undergo repeated asexual cell division in the host erythrocytes, during which parasite-coding proteins do extensive host cell modification. 2 The parasite extensively modifies the host cell membrane, resulting in changes in morphology, deformability, and adhesive properties of the host erythrocyte during pathogenesis. 3 Post malarial infection, the host erythrocytes become hard, reflected by the changes in the host membrane cytoskeleton structure. In addition to host cell modifications, parasite cytoskeletal machinery contributes to various biological processes such as parasite invasion, egression, cell motility, cell division, vesicle trafficking, mitochondrial morphology, cell polarity establishment, and gene regulation. 4, 5 Thus, the proteins involved in the cytoskeleton machinery can help better understand the related pathways and provide insights into the molecular details of important events accompanying parasite invasion to RBCs, cell division, and egression.

Extensive experimental evidence showed that microtubules associated with cytoskeleton machinery consist of the most tangled proteins involved in the pathogen's structural assembly, such as actin(s), tubulin(s), and histone(s). 6−8 P. falciparum consists of five tubulin superfamily members. The first two, namely, alpha (α)-tubulin and beta (β)-tubulin polymerize to form microtubules. Gamma (γ)-tubulin nucleates the microtubules as a component of the gamma-tubulin ring complex and is conserved. The delta (δ)-tubulin and epsilon (ε)-tubulin are majorly associated with centrioles, expressed in gametocytes, and their molecular function is unclear to date. 9 Among all the tubulins, Plasmodium encodes two isomers of αtubulin and a single copy of β-, γ-, δ-, and ε-tubulins, of which α2 and β are essential for parasite survival. 8−12 Genetic studies have also revealed that the presence of tubulins is essential for functions, which include microtubule nucleation, assembly, and polarity establishment. 8, 13 The basic structure of microtubules is the combination of αand β-tubulins that form heterodimers that bind head to tail during the protofilament development ( Figure 1A ). The γ-tubulin is important in the nucleation and polar orientation of microtubules. Microtubules also control the process of assembly and disassembly through guanine nucleotides bound to α1-, α2-, β-, and γ-tubulins. 14 The guanine triphosphate (GTP) and guanine diphosphate (GDP) are bound to these tubulins, which favors microtubule growth, and its hydrolysis results in shrinkage ( Figure 1B) . Tubulins have been shown to bind two GTP molecules, one at the exchangeable site of beta chains and another at a nonexchangeable site on the alpha chain. The P. falciparum α1tubulin (Pf3D7_0903700), α2-tubulin (Pf3D7_0422300), βtubulin (Pf3D7_1008700), and γ-tubulin (Pf3D7_0803700) are transcribed at all the stages of the parasite life cycle. Meanwhile, δ-tubulin (Pf3D7_0933800) and ε-tubulin (Pf3D7_1475700) are predominantly transcribed in gametocytes. 15 Although there is not much variation in amino acid sequences of α1-, α2-, β-, and γ-tubulins, however, there are vast differences in copy number variations, isotypes, and posttranslational modifications. 14,16 Therefore, focusing on the essentiality and microtubule polarity, the four tubulins are selected and marked as important targets in the current study. These four tubulins can be explored for novel drug discovery for malaria.

A recent study has identified a parasite-specific binding site in Plasmodium α-tubulin, similar to plant α-tubulin but absent in human α-tubulin. 17, 18 Recent advancements in structural bioinformatics, especially the methods based on comparative modeling, active site analysis, and protein simulations, have significantly progressed, which even facilitate the discovery of sites specific to a group of proteins. 19 A potential target must be essential for the pathogen, and further designing of inhibitors should hinder the function exclusively to pathogen and should avoid the undesirable cross-reactivity with the human proteins. The availability of complete genome sequences of pathogens in combination with Bioinformatics tools and databases is one of the great assistances to identify potential drug targets/vaccine candidates in a large pool of gene/protein polls. In this regard, based on essentiality, the four tubulins were marked as potential drug targets in microtubule machinery. The current study aims to design quality 3D structures of P. falciparum α1-, α2-, β-, and γtubulins by various homology modeling algorithms and ab initio methods. 20, 21 Additionally, molecular dynamics (MD) simulations were performed and atomic trajectories were analyzed by root-mean-square deviation (RMSD), root-meansquare fluctuation (RMSF), radius of gyration (Rg), solventaccessible surface area (SASA), secondary structure, principal component analysis (PCA), free-energy landscape (FEL), and density distribution analysis to evaluate the stability of the modeled tubulin structures. 22, 23 A common parasite-specific tubulin binding site embedded with GTP and GDP was also identified, which could be exploited for new antimalarial drug discovery targeted to the tubulin machinery. The analysis revealed that Leu92 and Arg228 are parasite-specific, and the amino acid residues Ala12, Ile14, Gln91, Ala100, Tyr179, and Ala180 present in α1-, α2-, and β-tubulins are parasite-specific too. Summarily, we successfully performed atomic resolution homology modeling, identifying a common parasite-specific interaction cleft, and MD simulation, which can potentially aid to identify novel antimalarial drugs, inhibitors, and vaccine candidates that can act toward the four tubulins' binding pocket. We believe that our study provides novel insights into exploring tubulins as novel next-generation therapeutic targets for malaria.

Conserveness of α1-, α2-, β-, and γ-Tubulins. The FASTA sequences of all the four P. falciparum tubulin proteins α1: Q6ZLZ9, α2: Q8IFP3, β: A0A143ZWL7, and γ: Q8IAN7 were retrieved from the UniProt database. The ClustalOmega was used for MSA; results revealed the highest scoring match and sequence identity between α1 (452), α2 (445), β (453), and γ (450) amino acids ( Figure 2 ). As expected, it shows a high similarity between the α1-, α2-, β-, and γ-tubulin sequences as they arise from the same family during evolution. The conserveness and the adaptive evolutionary patterns also strongly indicate the importance of the proteins during disease evolution ( Figures S1 and S2) .

Structure Prediction Outcomes. MODELLER. We constructed 3D structures for the α1-, α2-, β-, and γ-tubulins of P. falciparum to model the missing amino acids in the templates using MODELLER. The reliable templates were selected manually with BLASTp with identity ≥50% and query coverage ≥90% and aligned using CLUSTALX. After the successful run, 20 best models were generated based on DOPE and GA341 scores. The best models for all the α1-, α2-, β-, and γ-tubulins with the highest sequence identity score were selected as reliable models for respective proteins. The models B99990019 for α1-tubulin, B99990020 for α2-tubulin, B99990002 for β-tubulin, and B99990009 for γ-tubulin with the highest negative DOPE scores of −46193.89, −53800.18, −51552.94, and −26685.67 kcal/mol, respectively, were chosen for further analysis because of the measured conformational energy and relative structural stability ( Table 1 ). The details of selected templates, DOPE scores, identity, and query coverage generated using SWISS-MODEL, I-TASSER, and PHYRE2 are given in Table S1 .

SWISS-MODEL. The SWISS-MODEL generated the best alignments between the target-template sequences and homology models for α1-, α2-, β-, and γ-tubulins. Suitable templates with an identity of >60% were ranked as top 10, and the best models with higher GMQE and QMEAN scores were obtained with acceptable alignment values. The details of the templates selected, QSQE, GMQE, and QMEAN scores, along with identity and coverage, are shown in Table S2 . Moreover, evaluation of PROCHECK analysis for α1-, α2-, β-, and γtubulins revealed the evidence of their acceptability within the nominal range ( Figure S3 ).

I-TASSER. The FASTA sequences of the α1, α2, β, and γ proteins were used as inputs to the I-TASSER server, and the best template was chosen. The server analyzed the models for unstructured regions that correspond to the target protein's disordered areas using the ab initio modeling technique. The confidence scores (C-scores) of modeled structures of α1-, α2-, β-, and γ-tubulins were noted to be 1.04, 1.14, 1.10, and 1.28, respectively. The C-scores, cluster size, and cluster density are in the acceptable range, reflecting the good quality of the models. Further, the RMSD (1.0 Å) of the best first model for α1-, α2-, β-, and γ-tubulins was the least, and hence the models were selected for further analysis ( Figure S4 ). PHYRE2. The query sequence of the target proteins was screened against the curated sequences from PDB with >50% identity with the inbuilt database HHblits. PSI-BLAST-based secondary structure prediction was applied, and MSA was performed for α1-, α2-, β-, and γ-tubulin targets, which were combined into a query Hidden Markov Model. The results revealed top 10 models with high confidence scores. The top models for α1-, α2-, β-, and γ-tubulins were selected ( Figure  S5 ).

Validation of Predicted Homology Models. The validation of MODELLER-, SWISS-MODEL-, I-TASSER-, and PHYRE2-generated homology models was performed to identify the stereochemical properties of the residues in allowed regions for α1-, α2-, β-, and γ-tubulin. The results revealed that 90−98% of the residues are within the most favored regions, 8−10% residues of modeled proteins are within the additional allowed regions, 4−5% residues of modeled proteins are within the generously allowed regions, and only 0−2% residues were seen in disallowed regions. The stereochemical property analyses of the generated models reveal that the models are highly reliable in the distribution of φ and ψ angles. The predicted models of SWISS-MODEL, I-TASSER, and PHYRE2 programs (Table S2 ) revealed that the model generated with MODELLER is more acceptable due to the structures' stability via distribution of φ and ψ angles ( Table 1) .

Ramachandran Plot Analysis (RCP). The validations of the four reliable final models were achieved using PROCHECK, RAMPAGE, and MolProbity. The 3D structures constructed using the MODELLER algorithm for α1-tubulin revealed that 96% of residues were in allowed, generously allowed, and favored regions ( Table 1 ). The results also revealed that 0.8% of residues were in disallowed regions and only three residues were out of the box ( Figure 3A (i)). The results assigned by the RCP analysis for α1-tubulin obtained through RAMPAGE and MolProbity represented the best quality ( Figure 3B (i)). The generated structure of α2-tubulin protein revealed that 99% of residues were in the allowed region and favored regions, and only 0.5% of residues were in disallowed regions, representing the reliable and best-predicted homology model ( Figure  3A (ii),B(ii)). The 3D structure of β-tubulin using MOD-ELLER showed that 100% of residues were in allowed, generously allowed, and favored regions ( Figure 3A (iii),B(iii)). The plot statistics revealed that the β-tubulin model was highly reliable in quality. The results show that 98% of residues were in allowed regions and 1.2% of residues were in disallowed regions, representing the reliability of the model ( Figure  3A (iv),B(iv)). The RCP analysis performed using RAMPAGE and MolProbity revealed that the constructed γ-tubulin model has good stereochemical distribution. The detailed assessment of RCP analysis constructed using PROCHECK for the α1-, α2-, β-, and γ-tubulin models obtained by SWISS-MODEL, I-TASSER, and PHYRE2 programs is shown in Table S2 .

ProQ and ProSA Analyses. To evaluate the overall stereochemical properties and correctness of the four tubulin models, we used PROCHECK, ProQ, ProSA, and QMEAN analyses. Amongst the four modeling algorithms, the MODELLER results showed the best distribution of φ and ψ angles for α1-, α2-, β-, and γ-tubulins. The Z-scores by ProSA-predicted energy-minimized PDB structure values were −9.69, −9.72, −9.72, and −10.06 for the models for α1-, α2-, β-, and γ-tubulins, respectively ( Table 1 ). The Z-scores by SWISS-MODEL, I-TASSER, and PHYRE2 were also in the acceptable range, as shown in Table S3 . ProQ was applied and obtained reliable results with LGscores and MaxSub scores to understand structural features based on the neural network for model quality. The LGscores of α1-, α2-, β-, and γ-tubulins implied −log of P-values 7.346, 7.097, 6.423, and 6.518, respectively. Similarly, MaxSub values obtained were 0.673, 0.695, 0.575, and 0.568 for α1-, α2-, β-, and γ-tubulins, respectively, by ProQ, which show that all the values are in a significant range and reflect the reliability of the models ( Table  1 ). The LGscore and MaxSub score values for all the models generated by SWISS-MODEL, I-TASSER, and PHYRE2 algorithms are in an acceptable range (Table S3 and Figure  S6 ).

Binding Pocket Detection. The result of PDB-Sum, LigPlot, FunFold, 3D ligand site, Pocket-identifier, DoGSi-teScorer, and CastP tools found a common active pocket bound with GTP and GDP for α1, α2, β and γ modeled tubulins ( Figure S7 ). The α1-tubulin active pocket firmly embedded with GTP/GDP and Gln11, Ala12, Ile14, Gln91, Leu92, Ala100, Gly143, Gly144, Thr145, Tyr179, Ala180, and Arg228 amino acid residues. In the active cleft α2-tubulin bound with GDP/GTP, the residues, namely, Gln11, Ala12, Ile14, Gln91, Leu92, Gly143, Gly144, Thr145, and Arg228 are seen. The β-tubulin binding cavity site residues were Ala10, Gly11, Ile14, Gln91, Leu92, Ala12, Ala97, Asn99, Gly141, Gly142, Thr143, Gly144, and Arg228 along with GTP/GDP. Similarly, γ-tubulin binding site residues were Gly11, Ile14, Gln91, Leu92, Ala100, Asn102, Gly143, Thr145, Gly146, Pro173, Tyr179, Ala180, and Arg228 with the presence of GDP/GTP. From the interactive site analysis, it was evident that the α1-, α2-, β-, and γ-tubulins of P. falciparum were sharing a common binding pocket with amino acid residues Gln11, Ala12, Ile14, Gln91, Leu92, Ala100, Thr143, Gly144, Tyr179, Ala180, and Arg228 with the presence of GTP and GDP (Figures S8−S10). A recent study has shown that residues Arg2, Gln133, Arg243, Asn249, Val250, Asp251, Val252, Thr253, and Glu254 are parasite α-tubulin-specific, 17 and a detailed comparative analysis of the four tubulin proteins with human tubulins revealed that Leu92 and Arg228 are found to be parasite-specific, which were lying around the same active region. Additionally, we found that the amino acid residues Ala12, Ile14, Gln91, Ala100, Tyr179, and Ala180 present in α1-, α2-, and β-tubulins are parasite-specific too. The diagrammatic representation of the identified parasitespecific residues and the interaction pocket of the amino acid residues is depicted in Figure 2 , Figure 4 , and Table S4 .

MD Simulation Interpretation. We used GROMACS in a Linux environment to run MD simulations. We analyzed the obtained MD trajectories to study parameters like RMSD, RMSF, Rg, SASA, intramolecular hydrogen bond, PCA, FEL, and density distribution analysis. 23, 24 The MD simulation studies were carried out to understand α1-, α2-, β-, and γtubulin dynamic behaviors in the presence of water and examine the protein stability.

RMS Deviations and RMS Fluctuations. The deviations were calculated on four tubulins defined from the starting point of simulation to the end frame via a consistent map of RMS deviation and fluctuation plots over the time scale. We found that the average RMSD value of α1-tubulin was 0.648 ± 0.01 nm, α2-tubulin was 0.277 ± 0.01 nm, β-tubulin was 0.588 ± 0.01 nm, and γ-tubulin was 0.430 ± 0.01 nm ( Figure 5A (i− iv)). The plateau of the average RMSD values was directly proportional to the modeled tubulins' average stable conformations (Table S5 ). The graph portrayed that the α1tubulin reached equilibrium from 20 ns, α2-tubulin reached a steady-state from 30 ns, β-tubulin displayed stability from 10 ns, and γ-tubulin presented a steady state from 30 ns. Moreover, the initial and the time-framed RMSD in the last trajectory were noticeably stable for the α1-, α2-, β-, and γtubulin models.

The RMS fluctuations for the four tubulins from their timeaveraged values are stable, and there are no major differences between the start and end time of simulation trajectories. The RMSF analysis revealed that the average values of α1-, α2-, β-, and γ-tubulins were within the nominal ranges by sharing common binding pocket residues, namely, Leu92, Arg228, Ala12, Ile14, Gln91, Ala100, Tyr179, and Ala180, corroborating with the interaction site analysis ( Figure 5B ). The initial and the end fluctuations at some residues may be discounted owing to relaxation in the homology model induced by the force field. From the analysis, it is evident that the average distances between the atoms from the modeled α1-, α2-, β-, and γ-tubulin structures demonstrate nominal residue displacement within the flexible range.

Determination of Rg, SASA, and Hydrogen Bond Analyses. The Rg analysis of the four tubulins defined the mass-weighted root-mean-square distance of atoms from their respective center of mass. The calculated average Rg was found to be nominal, with the values of α1-tubulin = 2.29 ± 0.01 nm, α2-tubulin = 2.16 ± 0.01 nm, β-tubulin = 2.09 ± 0.01 nm, and γ-tubulin = 2.74 ± 0.01 nm ( Figure 6A (i−iv) and Table S5 ). Moreover, the overall dimensional stable compactness of the modeled structures had high compactness between the residues and the bonding pattern, in turn portraying a direct proportion to the protein volume.

The exposed surface area of modeled tubulin structures was subjected to SASA analysis. The SASA clearly illustrates that the protein solvent-accessible surface area is directly proportional to the tendency of the interaction of self-fixed and other molecules. The rearrangement of the hydrogen bond network between the side chain residues and the surrounding water molecules of the tubulins' modeled structure was calculated. The average SASA of all the protein was noticeably different in values, with 274.43 ± 0.01 nm for α1-tubulin, 262.89 ± 0.01 nm for α2-tubulin, 242.87 ± 0.01 nm for β-tubulin, and 274.90 ± 0.01 nm for γ-tubulin ( Figure 6B(i−iv) and Table S5 ). The SASA results depicted that the solvent surface area of the four proteins was highly flexible.

The stability of the protein is determined by the number of hydrogen bonds formed between the atoms. The numbers of H-bonds were calculated to understand the stability of the four tubulins. The hydrogen bond analysis revealed the average hydrogen bonds to be 281.15, 321.75, 352.69, and 310.52 bonds for α1-, α2-, β-, and γ-tubulins, respectively ( Figure  6C(i−iv) ). The four tubulins displayed a high number of hydrogen bonds, with no alterations in forming the intramolecular bond numbers, revealing strong contacts between the residues in the protein.

Secondary Structures and Principal Component Analysis (PCA). The secondary structures of the tubulins were investigated to check the stability in the formation of helix, sheets, coils, loops, and turns to understand the structural plasticity. We found that the common active site residues of all tubulins fall in alpha-helix regions, which include Ala12, Ile14, Gln91, Leu92, Ala100, Tyr179, Ala180, and Arg228. The secondary structure analysis demonstrated that the percentages of the structural elements in the modeled structures were constant for the tubulins over time at each nanosecond (ns) ( Figure 7A(i−iv) ).

The conformational subspace of the complexes is examined to understand the four tubulin models' stable dynamic behavior by the PCA technique. The results revealed that the clusters are well defined and covered minimum subspaces on the proteins. The graph was plotted between the eigenvectors 1 and 2, which gave the displacement of atomic fluctuation and motions in the structures. The values ranged from −4 to 4 nm 2 for eigenvector 1 and −6 to 7 nm 2 for eigenvector 2 for the α1-tubulin model. The α2-tubulin showed eigenvector 1 values of −4 to 4 nm 2 and eigenvector 2 values of −6 to 5 nm 2 of minimum cluster range. β-tubulin revealed eigenvector 1 values of −4 to 6 nm 2 and eigenvector 2 values of −6 to 10 nm 2 , and γ-tubulin occupied a small occupancy cluster range with eigenvector 1 values of −4 to 3 nm 2 and eigenvector 2 values of −6 to 2 nm 2 , respectively ( Figure 7B(i−iv) ). The eigenvector ranges of all the four tubulin models showed restricted space, leading to a welldefined internal motion behavior, vital for model stabilization.

Free-Energy Landscape (FEL) and Density Distribution Analyses. We performed FEL analysis of the modeled tubulin, which accurately describes the protein folding behavior. 23, 25 The protein attaining a specific position can be marked by the consumption of total energy on the landscapes. The global free-energy minima were 11.7 kJ/mol for α1-tubulin, 13.8 kJ/ mol for α2-tubulin, 14.4 kJ/mol for β-tubulin, and 13.5 kJ/mol for γ-tubulin ( Figure 8A(i−iv) ). The values of multiple energy minima revealed stable significance for the four modeled structures of tubulins, revealing a stable conformation. The analysis also confirmed to have a stable folding pattern in all the homology modeled α1-, α2-, β-, and γ-tubulin structures.

To understand the atomic density, atomic orientation, and distribution of modeled tubulin structures, we performed density distribution analysis of the molecular coordinates of each of the models during MD simulations. 23, 24 The stable density area was observed for α1-tubulin with a value of 3.49 nm −3 , α2-tubulin with a value of 5.27 nm −3 , β-tubulin with a value of 3.87 nm −3 , and γ-tubulin with a value of 4.64 nm −3 ( Figure 8B(i−iv) ). The density area analysis confirms that the tubulin structures are stable with minimum energy wells.

Evidence from genetic studies of Plasmodium has revealed that the presence of tubulins is essential in microtubule nucleation, assembly, and polarity establishment. All the α-, β-, and γtubulins are considered vital targets due to their essentiality and presence in all the parasite life cycle stages. In this regard, we have constructed different 3D homology models of P. falciparum α1-, α2-, β-, and γ-tubulins. The study enabled the generation of energy-refined and reliable quality 3D models based on homologous tubulin structures based on protein homology modeling, stereochemical assessment, and thermo- dynamic evaluations. MODELLER-generated models displayed the best distribution of φ and ψ angles for α1-, α2-, β-, and γ-tubulin models.

Consequently, we performed MD simulations for the tubulins to understand the stability and dynamics of the models. The RMSD, RMSF, Rg, SASA, PCA, FEL, and density distribution analysis reciprocated with the simulation behavior. We also determined a common binding cleft and parasitespecific amino acid residues in the four tubulins when compared to human sequences, which can facilitate the designing of novel inhibitors specifically targeting the malaria parasite-specific tubulin machinery.

Hardware and Software. The work was carried out on a High-Performance Computing (HPC) IBM server with 128GB RAM in each node and four CUDA-enabled NVIDIA (model: Nvidia Tesla V100) graphics processing units (GPUs) with 16GB RAM each. In the current study, we utilize various online tools such as CASTp, PlasmoDB, PHYRE2, PRO-CHECK, RAMPAGE, MolProbity, FunFold, DoGSiteScorer, Pocket-identifier, PDB-Sum, LigPlot, ProQ, ProSA, ClustalOmega, SWISS-MODEL, and I-TASSER.

Selection of Tubulins. PlasmoDB is a widely used database by malaria researchers that consists of the emerging completed genome sequences and annotation of sequencing projects on Plasmodium spp. 15, 26 The sequence information of proteins is integrated with other genomic-scale data of the Plasmodium research community and includes gene expression data, microarray projects, and proteomics studies.

Multiple Sequence Alignment (MSA). MSA is often used to assess conserved sequences of protein domains, secondary and tertiary structures, and even individual amino acids or nucleotides. 27,28 CLUSTAL (X-online and W-offline) is a widely used MSA computer program. The tubulin sequences were aligned using ClustalOmega 29 and CLUS-TALX 30 to prepare the initial alignment files. CLUSTALX, the CLUSTAL software with an offline graphical user interface, was utilized to align the sequences and study their conservation among the tubulins. The tubulin sequences were subjected to ClustalOmega to understand the percentage conservation, and the results were cross-validated with CLUSTALX. 30 Target-Template Alignment. The FASTA formatted sequences of the selected tubulins were retrieved from the UniProt database by providing a 6-digit unique ID as inputs. The target sequences were queried using Basic Local Alignment Search Tool protein (BLASTp, https://blast.ncbi. nlm.nih.gov/) against the Protein Data Bank (PDB) sequences 31 to select the templates. The proteins with high identity, similarity, and query coverage of ≥50%, ≥70%, and ≥90% were selected as templates, respectively. The best selected target-template sequences were aligned manually using MSA and subjected to protein modeling algorithms to construct a reliable 3D model.

Comparative Homology Modeling Algorithms. In the present study, MODELLER 32 was used to perform homology modeling for the selected tubulin protein sequences. We also performed comparative homology modeling using various modeling algorithms such as I-TASSER, 33 PHYRE2, 34 and SWISS-MODEL. 35 The salient details of each of the modeling algorithms is given below:

MODELLER. The software package is the most popular and widely used molecular modeling technique to predict protein 3D structures on the basis of sequence similarity. The package provides model evaluation to minimize modeling errors. In the case of more than one model, the best model is chosen using Discrete Optimized Protein Energy (DOPE) assessment and GA341 score (displayed at the terminal and in the log file) at the end of runs. The MODELLER scripts for generation of models were developed based on heteroatom modeling with 20 models. The best 3D structures were selected for further validation.

SWISS-MODEL. The SWISS-MODEL workspace is freely accessible at http://swissmodel.expasy.org/workspace/. The algorithm is the first publicly available and widely used automated modeling server. The amino acid sequence in FASTA format of the target protein is used as inputs. The algorithms perform BLASTp and provide the best homology models based on Global Model Quality Estimation (GMQE) and Qualitative Model Energy ANalysis (QMEAN) statistical parameters. The GMQE provides the quality by combining the properties from the target-template alignment and estimates a range between 0 and 1. The scoring function of QMEAN consists of a linear combination of structural descriptors, which is related to the high-resolution X-ray crystal structure and estimate Z-score scheme.

Iterative Threading ASSEmbly Refinement (I-TASSER). The technique combines different approaches with a framework of reassembly and ab initio methods. The Monte Carlo algorithm is implemented in the I-TASSER suite, where modeling is performed for the initial models by a knowledge-based energy function with spatial constraints and hydrogen bonds. 36 The structural templates are identified first from PDB through multiple threading and include full-length atomic models generated with iterative template-based fragment assembly simulations. The best quality model from the I-TASSER provides good confidence scores and RMSD values.

Protein Homology/AnalogY Recognition Engine (PHYRE2). This method is a spontaneous and straightforward interface web tool to predict the protein structure analysis using the threading technique (http://www.sbg.bio.ic.ac.uk/ phyre2). Apart from building the protein model, it can predict the ligand-binding sites and analyze amino acid variants for a given protein query. PHYRE2 can also determine the secondary structures, domain composition, and quality of the model.

Validation of the Best Models. The 3D model validation of the tubulin models was performed using a Ramachandran plot (PROCHECK), 37 Protein Quality predictor (ProQ), 38 RAMPAGE (http://mordred.bioc.cam.ac.uk/~rapper/ rampage.php), MolProbity, 39 and Protein Structure Analysis (ProSA). 40 The Ramachandran plot of 3D structures can be determined using torsion angles for every individual amino acid of the protein. The torsion angle should fall in the allowed and favored regions with more than 90% occupancy for the best quality model. The plot statistics obtained from RCP reveals that when the residues are within allowed areas with high resolution, then the structures tend to have good clustering. The best 3D model was verified using ProQ and ProSA to check the generated models' potential errors. Both the web servers validate the overall quality score depending on the comparison of input structures with high-resolution experimentally refined protein structures with local quality estimation, Z-score, LGscore, and MaxSub cutoff. For a reliable and extreme quality model with Z-scores > 4.0, LGscore >3 is good and >5 is a very good model and MaxSub >0.5 is good and >0.8 is a very good model.

Identification of Interaction Pockets. The binding pocket for protein/enzyme is important to understand the molecular and biological functions. To analyze the favorable binding between the protein and ligand, various bioinformatics prediction tools were outstretched. The binding cleft and the interaction site amino acid analysis of the tubulins were predicted using FunFold, 41 3D ligand site, PDB-Sum, 42 DoGSiteScorer, 43 and CASTp tools. 44 The interactions between protein and the cocrystal molecules were plotted using LigPlot 45 and visualized using PyMol. 46 MD Simulations for the Modeled Structures. GROningen MAchine for Chemical Simulations (GROMACS) 47 is a package to perform and analyze MD simulations. MD simulations were carried out using a GROMACS 5.18.3 software package. 47−49 The modeled 3D structures of tubulin were utilized for the MD simulations. The α1-, α2,-β-, and γtubulins' topology parameter files were created by the CHARMM27 force field included with CMAP correction. 50 The intermolecular (nonbonded) potential, represented as the sum of Lennard-Jones (LJ) force-based switching with a cutoff distance range of 8−10 Å, pairwise Coulomb interaction, and long-range electrostatic force were determined by a particle mesh Ewald (PME) approach. 51 The real-space cutoff was set to 1.2 nm for the PME calculations too. 52 The velocity Verlet algorithm was applied for the numerical integrations and the initial atomic velocities were generated with a Maxwellian distribution at the given absolute temperature. The system was then immersed with the default TIP3P water model, and protein was placed at the center of the cubic grid box (1.0 nm 3 ). 53 Neutralization was performed to make the concentration of the system 0.15 M. The neutralized system was then subjected to energy minimization using the Steepest Descent (SD) and Conjugate Gradient (CG) algorithms utilizing a convergence criterion of 0.005 kcal mol −1 Å −1 . The twostandard equilibration phase was carried out separately under NVT (atom, volume, temperature) and NPT (atom, pressure, temperature) ensemble conditions such as constant volume and constant pressure for each complex similar time scale with LINear Constraint Solver (LINCS) to restrain the bonds involving hydrogen atoms. A Nose−Hoover thermostat and Parinello−Rahman barostat were applied to maintain the temperature and the pressure of the system. The system was maintained constant at 1 bar and 300 K, with a coupling time of τP = 2 ps and τT = .1 ps. The Periodic Boundary Condition (PBC) is used for integrating the equation of motion by applying the leap-frog algorithm with a 2 fs time step. Finally, to make the system ready for production, the fixing of constraints is achieved with the relaxation of the grid box with water and counterions. 54, 55 The simulation was performed for 100 ns run for the modeled α1-, α2-, β-, and γ-tubulin structures and all the trajectories were recorded.

• The Cα-atom fluctuations were calculated by RMSD with the equation: where RMSF i is the RMSF of i atom and r i is the atom position of i. The lower the RMSF, the higher the rigidity, while the higher the RMSF, the lower the flexibility.

• The Rg was calculated to measure the compactness, which is given by the equation The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acsomega.1c01988. Phylogeny evolutionary relationship between the four tubulin proteins; conserveness analysis of the four tubulin proteins of P. falciparum and templates; Ramachandran plot analysis of α1-, α2-, β-, and γtubulin targets using the SWISS-MODEL prediction algorithm; Ramachandran plots obtained through PDB-Sum analysis for the tubulin targets using I-TASSER modeling algorithms; PROCHECK-obtained stereochemical analysis for the selected four tubulin targets using the PHYRE2 modeling prediction method; quality check using ProSA showing depiction of Z-scores and local quality models for the selected tubulin proteins; homology modeled reliable structures of the respective four tubulin targets from the MODELLER prediction method visualized through PyMOL; active site amino acid display of the selected four tubulin targets obtained through LigPlot analysis; interactive site representation of amino acids between GTD, GDP, and the four tubulin proteins obtained through Fun-Fold prediction server; representation of α1-, α2-, β-, and γ-tubulin active sites identified around GTD and GDP molecules aligned with overlapping module of PyMOL; templates selected with different modeling algorithms; Ramachandran plot analysis with other modeling tools; Z-scores,

LGscores, and MaxSub scores of predicted models of tubulin by ProSA and ProQ quality checking tools; interactive amino acid site of tubulin targets embedded with GTP and GDP molecules; and average stability scores of the predicted four tubulin targets obtained during the 100 ns simulation run time (PDF) 

World malaria report 2020: 20 years of global progress and challenges

Plasmodium asexual growth and sexual development in the haematopoietic niche of the host

Identification of essential exported Plasmodium falciparum protein kinases in malariainfected red blood cells

In silico designing and analysis of inhibitors against target protein identified through host-pathogen protein interactions in malaria

The parasitophorous vacuole of the blood-stage malaria parasite

Functional homo-and heterodimeric actin capping proteins from the malaria parasite

A histone methyltransferase inhibitor can reverse epigenetically acquired drug resistance in the malaria parasite Plasmodium falciparum

Functional profiling of a Plasmodium genome reveals an abundance of essential genes

Interaction of Plasmodium falciparum apicortin with α-and βtubulin is critical for parasite growth and survival

α-Tubulin II is a male-specific protein in Plasmodium falciparum

The γ-tubulin gene of the malaria parasite Plasmodium falciparum

A protein interaction network of the malaria parasite Plasmodium falciparum

Zeta-tubulin is a member of a conserved tubulin module and is a component of the centriolar basal foot in multiciliated cells

New insights into microtubule elongation mechanisms

PlasmoDB: exploring genomics and post-genomics data of the malaria parasite, Plasmodium falciparum

An epigenetic map of malaria parasite development from host to vector

Alanine scanning of dinitroaniline/ phosphorothioamidate site of α-tubulin in plasmodium species distributed in India

In-depth comparative analysis of malaria parasite genomes reveals protein-coding genes linked to human disease in Plasmodium falciparum genome

An overview of comparative modelling and resources dedicated to large-scale modelling of genome sequences

SWISS-MODEL: an automated protein homology-modeling server

Structural and functional alterations of nitric oxide synthase 3 due to missense variants associate with high-altitude pulmonary edema through dynamic study

Screening Malaria-box compounds to identify potential inhibitors against SARS-CoV-2 Mpro, using molecular docking and dynamics simulation studies

Insights into the structural and dynamical changes of spike glycoprotein mutations associated with SARS-CoV-2 host receptor binding

Simultaneous Inhibition of SARS-CoV-2 Entry Pathways by Cyclosporine

Molecular motions and free-energy landscape of serine proteinase K in relation to its cold-adaptation: a comparative molecular dynamics simulation study and the underlying mechanisms

PlasmoDB: the Plasmodium genome resource. A database integrating experimental and computational data

Multiple sequence alignment with Clustal X

202 Subunit vaccine design against pathogens causing atherosclerosis

Clustal Omega, accurate alignment of very large numbers of sequences

Clustal W and Clustal X version 2.0

The protein data bank

Evaluation of comparative protein modeling by MODELLER

I-TASSER server for protein 3D structure prediction

The Phyre2 web portal for protein modeling, prediction and analysis

SWISS-MODEL and the Swiss-Pdb Viewer: an environment for comparative protein modeling

The I-TASSER Suite: protein structure and function prediction

Structure validation by Cα geometry: ϕ, ψ and Cβ deviation

Can correct protein models be identified?

MolProbity: More and better reference data for improved all-atom structure validation

ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins

The FunFOLD2 server for the prediction of protein−ligand interactions

PDBsum: summaries and analyses of PDB structures

DoGSiteScorer: a web server for automatic binding site prediction, analysis and druggability assessment

CASTp 3.0: computed atlas of surface topography of proteins

LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions

High performance molecular simulations through multi-level parallelism from laptops to supercomputers

5: a high-throughput and highly parallel open source molecular simulation toolkit

Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation

CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields

OpenMM simulations using the CHARMM36 additive force field

A critical appraisal of the zero-multipole method: Structural, thermodynamic, dielectric, and dynamical properties of a water system

III A modified TIP3P water potential for simulation with Ewald summation

Designing of phenolbased β− carbonic anhydrase1 inhibitors through QSAR, molecular docking, and MD simulation approach

2/3D-QSAR, molecular docking and MD simulation studies of FtsZ protein targeting benzimidazoles derivatives