Atomic Resolution Homology Models and Molecular Dynamics Simulations of Plasmodium falciparum Tubulins.

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.

Introduction

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 A). 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 B). Tubulins have been shown to bind two GTP molecules, one at the exchangeable site of beta chains and another at a non-exchangeable site on the alpha chain. The P. falciparum α1-tubulin (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 post-translational 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.

Figure 1

Schematic diagram of tubulin machinery. (A) Heterodimer complex formation of α- and β-tubulins. (B) Representation of assembly and disassembly processes by α-, β-, and γ-tubulins with GTP and GDP involvement.

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-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible 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.

Results and Discussion

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 ). 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).

Figure 2

MSA of α1-, α2-, β-, and γ-tubulins and templates along with PlasmoDB, UniProt, and PDB IDs and the common binding pocket and the parasite-specific amino acid residues (92 and 228) are highlighted. (α1) target: Pf3D7_0903700, template: 5JC0; (α2) target: Pf3D7_0422300, template: 6E7B; (β) target: Pf3D7_1008700, template: 6NNG, and (γ) target: Pf3D7_0803700, template: 3CB2.

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 ). The details of selected templates, DOPE scores, identity, and query coverage generated using SWISS-MODEL, I-TASSER, and PHYRE2 are given in Table S1.

Table 1

Template Details, Ramachandran Plot Analysis, Z-Scores, LGscores, and MaxSub Scores of MODELLER Predicted Tubulin Models

Details of template
protein	template	identity	coverage	DOPE	GA341
α1-tubulin	5JCO	84.90%	96%	–46193.89 kcal/mol	1.00
α2-tubulin	6E7B	83.83%	97%	–53800.18 kcal/mol	1.00
β-tubulin	6NNG	95.35%	96%	–51552.94 kcal/mol	1.00
γ-tubulin	3CB2	64.52%	99%	–26685.67 kcal/mol	1.00

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 ).

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 ). The results also revealed that 0.8% of residues were in disallowed regions and only three residues were out of the box (Figure A(i)). The results assigned by the RCP analysis for α1-tubulin obtained through RAMPAGE and MolProbity represented the best quality (Figure B(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 A(ii),B(ii)). The 3D structure of β-tubulin using MODELLER showed that 100% of residues were in allowed, generously allowed, and favored regions (Figure A(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 A(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.

Figure 3

Stereochemical properties and distribution of φ and ψ angles in tubulins. (A) RAMPAGE and (B) MolProbity: (i) α1-tubulin, (ii) α2-tubulin, (iii) β-tubulin, and (iv) γ-tubulin.

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 ). 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 ). 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, DoGSiteScorer, 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 parasite-specific residues and the interaction pocket of the amino acid residues is depicted in Figure , Figure , and Table S4.

Figure 4

Representation of the common and active interaction pocket of respective α1-tubulin represented in salmon color, α2-tubulin in green, β-tubulin in cyan, and γ-tubulin in blue visualized by PyMOL.

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 A(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 α1-tubulin 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.

Figure 5

RMSD and RMSF plots. (A) Depiction of the RMSD plot at given time points having stable consistency in the modeled tubulin structures captured during protein simulations: (i) α1-tubulin, (ii) α2-tubulin, (iii) β-tubulin, and (iv) γ-tubulin at 100 ns each. (B) RMSF plot showing fluctuations at residual of similar binding pocket amino acids in all modeled tubulin structures.

The RMS fluctuations for the four tubulins from their time-averaged 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 B). 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 A(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.

Figure 6

Rg, SASA, and hydrogen bond during MD simulations. (A) (i–iv) Representation of Rg during 100 ns MD simulation for each of the modeled tubulin structures. (B) (i–iv) Analysis of SASA on the tubulin homology model during MD simulations. (C) (i–iv) Hydrogen bond numbers during MD simulations.

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 B(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 C(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 A(i–iv)).

Figure 7

Secondary structure elements and projection of eigenvectors 1 and 2. (A) (i–iv) Secondary structure elements of the generated tubulin modeled structures captured at the simulation time. (B) (i–iv) Dynamic energy fluctuation plotted between two projection on eigenvectors 1 and 2 generated for the modeled structures showing conformational space of Cα-atoms.

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 nm2 for eigenvector 1 and −6 to 7 nm2 for eigenvector 2 for the α1-tubulin model. The α2-tubulin showed eigenvector 1 values of −4 to 4 nm2 and eigenvector 2 values of −6 to 5 nm2 of minimum cluster range. β-tubulin revealed eigenvector 1 values of −4 to 6 nm2 and eigenvector 2 values of −6 to 10 nm2, and γ-tubulin occupied a small occupancy cluster range with eigenvector 1 values of −4 to 3 nm2 and eigenvector 2 values of −6 to 2 nm2, respectively (Figure B(i–iv)). The eigenvector ranges of all the four tubulin models showed restricted space, leading to a well-defined 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 A(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.

Figure 8

FEL and density distribution graphs. (A) The energy wells represent tubulin proteins to identify the folding behavior and direction of motions of modeled structures analyzed with FEL. (B) Depiction of the density distribution of tubulin structures. (i) α1-tubulin, (ii) α2-tubulin, (ii) β-tubulin, and (iv) γ-tubulin.

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 B(i–iv)). The density area analysis confirms that the tubulin structures are stable with minimum energy wells.

Conclusions

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 thermodynamic 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 parasite-specific 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.

Materials and Methods

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, PROCHECK, 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 CLUSTALX[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:

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.

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 nm3).[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 two-standard 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.M = Σ, and r(t) is the position of atom i at time t after least square fitting the structure to the reference structure.where RMSF is the RMSF of i atom and r is the atom position of i. The lower the RMSF, the higher the rigidity, while the higher the RMSF, the lower the flexibility.where m is the mass of atom i and r the position of atom i with respect to the center of the mass of the molecule.where *Nacc and Ntot are the number of the accessible and total number of points in the shell and R is the sum of van der Waals radius.

The Cα-atom fluctuations were calculated by RMSD with the equation:

The flexibility of the structures was characterized using the RMSF of their Cα-atoms, calculated by

The Rg was calculated to measure the compactness, which is given by the equation

SASA was defined by the following equation: