A Computational Study of Molecular Mechanism of Chloroquine Resistance by Chloroquine Resistance Transporter Protein of Plasmodium falciparum via Molecular Modeling and Molecular Simulations

The molecular mechanism of chloroquine resistance by the chloroquine resistance transporter protein of Plasmodium sp. is explored using molecular modeling and computational methods. The key mutation, lysine(K)-76 to threonine(T) (LYS76THR) in the transporter protein pertains to increased recognition of the protonated forms of the antimalarial drug. Such enhanced affinity can promote drug efflux from host digestive vacuole, rendering aminoquinoline-based treatment ineffective.


Introduction
Malaria, one of the deadliest parasitic diseases of tropical origin, is caused by the parasite from Plasmodium genus [1]. The most effective treatments of malaria involve administration of aminoquinoline-based anti-malarial drug(s) either as a standalone treatment or in combination therapy [2][3][4]. However, emergence of resistance against both aminoquinoline drugs as well as artemisinin in combination therapy, has rendered this treatment strategy severely handicapped [5][6][7][8]. The once wonder-drug chloroquine (CQ) is no longer recommended for the treatment of P. falciparum malaria. One of the highly cited mechanism of CQ involves binding to iron protoporphyrin-IX, which prevents heme detoxification, thus poisoning the parasite with its own metabolic waste. The mechanism of CQ resistance (CQR) is an open debate to date. CQR was primarily attributed by Fitch to the deficiency of CQ binding in the resistant P. berghei-infected cells [9]. The deficiency in CQ binding in CQR strains can either be attributed via involvement of a molecular machinery transporting protonated drug molecules out of the host digestive vacuole (DV) or via differential binding of different protonated forms of CQ to heme in CQ-sensitive (CQS) and CQR species [10]. CQ is a weak base and can freely diffuse across membrane in its neutral form and concentrate in the acidic DV. Inside DV, CQ protonates to membrane-impermeable bis-protonated CQH 2+ (Figure 1). Over the years, several genes associated with transporter proteins are attributed to quinoline resistance in P. falciparum. Important among these proteins are chloroquine resistance transporter (PfCRT), the multidrug resistance transporter 1 (PfMDR1), and the multidrug resistance-associated protein 1 (PfMRP1) [11][12][13][14][15]. While several studies have pointed out the multifaceted nature of the CQR, the most common resistance mechanism used by parasites to remove CQ is identified as mutation(s) in the chloroquine-resistant transporter of P. falciparum (PfCRT) protein and is thought to be involved in increased drug recognition and subsequent perturbation of drug concentration in the DV either via an efflux system [16,17] or by providing a path for protonated CQ to leak against concentration gradient [18]. Evidence of the ability of PfCRT to interact Physchem 2021, 1 233 with drugs was reported via fluorescence-based assay of the drug-associated efflux of protons from parasite DV [19,20]. Another proposal was that CRT is actually a chloride transporter, which altered food vacuole pH via chloride conductance [21,22]. There is a possibility that mutant PfCRT may affect the H + -transport across the DV membrane, potentially short-circuiting the H + -pump [23][24][25]. One interesting consequence of PfCRT mutation is swelling of DV, presumably due to osmotic pressure derived from inability of the mutant protein to transport its natural substrate(s) [26].
Physchem 2021, 1, FOR PEER REVIEW 2 fied as mutation(s) in the chloroquine-resistant transporter of P. falciparum (PfCRT) protein and is thought to be involved in increased drug recognition and subsequent perturbation of drug concentration in the DV either via an efflux system [16,17] or by providing a path for protonated CQ to leak against concentration gradient [18]. Evidence of the ability of PfCRT to interact with drugs was reported via fluorescence-based assay of the drugassociated efflux of protons from parasite DV [19,20]. Another proposal was that CRT is actually a chloride transporter, which altered food vacuole pH via chloride conductance [21,22]. There is a possibility that mutant PfCRT may affect the H + -transport across the DV membrane, potentially short-circuiting the H + -pump [23][24][25]. One interesting consequence of PfCRT mutation is swelling of DV, presumably due to osmotic pressure derived from inability of the mutant protein to transport its natural substrate(s) [26]. The identification of the CRT gene (MAL7P1.27) as the primary determinant of CQR was reported by Fidock et al. [21]. The pfcrt-gene encodes a 48-kDa putative transporter protein of 424 amino acid residues with high percentage of predicted transmembrane domains. Phylogenetic analysis predicts PfCRT to be a member of the drug/metabolite transporter superfamily of electrochemical potential driven transporters [27]. Additionally, PfCRT is an extraordinary polymorphic protein with no less than thirty variants known to date. In all resistant parasites, lysine-76 (76 K) is mutated to uncharged amino acids viz. threonine (76 T) or alanine (76 A) in field isolates [28], with the former shown to be the essential to show CQR [29]. Purification of PfCRT from yeast has been reported by two laboratories, although no structural information, either solution phase or solid state, is known [30,31]. The absence of a three-dimentional structure of PfCRT with a combination of a plethora of reports on its role in drug resistance has raised several speculations about the mechanism of CQ resistance. The structure and mechanism of chloroquine resistance is absent in the literature [32]. The present work aims to explore the underlying reasons for chloroquine resistance via mutation of PfCRT protein using computational models of wild type and mutant PfCRTs and on differential recognition of CQ.

Materials and Methods
Several bioinformatics analyses of structure and function of PfCRT place it either as a member of drug/metabolite transporter superfamily or as a gated aqueous pore [33]. Hence, we decided first to analyze the PfCRT sequence for structural similarity among known proteins. This exercise is aimed to find suitable template(s) to build homology models of PfCRT and further refining models thus obtained using molecular dynamics methods. Top models, wild type and mutated, were used for molecular docking calculations with both neutral and protonated CQ as ligands. The PfCRT-CQ complexes were subjected to MD simulations to explore the binding of the CQ ligand on mutant and wildtype variants of PfCRT.

Homology Model Building of PfCRT
The amino acid sequence of PfCRT is retrieved from UniPROT [34]. Our initial attempt to search for a suitable modeling template from a protein databank using the PSI- The identification of the CRT gene (MAL7P1.27) as the primary determinant of CQR was reported by Fidock et al. [21]. The pfcrt-gene encodes a 48-kDa putative transporter protein of 424 amino acid residues with high percentage of predicted transmembrane domains. Phylogenetic analysis predicts PfCRT to be a member of the drug/metabolite transporter superfamily of electrochemical potential driven transporters [27]. Additionally, PfCRT is an extraordinary polymorphic protein with no less than thirty variants known to date. In all resistant parasites, lysine-76 (76 K) is mutated to uncharged amino acids viz. threonine (76 T) or alanine (76 A) in field isolates [28], with the former shown to be the essential to show CQR [29]. Purification of PfCRT from yeast has been reported by two laboratories, although no structural information, either solution phase or solid state, is known [30,31]. The absence of a three-dimentional structure of PfCRT with a combination of a plethora of reports on its role in drug resistance has raised several speculations about the mechanism of CQ resistance. The structure and mechanism of chloroquine resistance is absent in the literature [32]. The present work aims to explore the underlying reasons for chloroquine resistance via mutation of PfCRT protein using computational models of wild type and mutant PfCRTs and on differential recognition of CQ.

Materials and Methods
Several bioinformatics analyses of structure and function of PfCRT place it either as a member of drug/metabolite transporter superfamily or as a gated aqueous pore [33]. Hence, we decided first to analyze the PfCRT sequence for structural similarity among known proteins. This exercise is aimed to find suitable template(s) to build homology models of PfCRT and further refining models thus obtained using molecular dynamics methods. Top models, wild type and mutated, were used for molecular docking calculations with both neutral and protonated CQ as ligands. The PfCRT-CQ complexes were subjected to MD simulations to explore the binding of the CQ ligand on mutant and wildtype variants of PfCRT.

Homology Model Building of PfCRT
The amino acid sequence of PfCRT is retrieved from UniPROT [34]. Our initial attempt to search for a suitable modeling template from a protein databank using the PSI-BLAST server at NCBI did not yield significant success. Subsequently, we used the PSIPRED server for fold recognition using the GenThreader method [35][36][37][38]. The templates with high and medium confidence, thus obtained, were used for the next step of homology modeling (vide supra). One interesting candidate is the glycerol phosphate transporter (PDB code: 1PW4), which is a membrane-bound transporter protein with a large area of Physchem 2021, 1 234 sequence match for homology modeling. A sequence alignment is then carried out using the Clustal-W2 program using the Blossum matrix. A total of 500 homology models were generated with the standard "automodel" class in MODELLER [39][40][41][42], and the model with the best DOPE score was selected for further refinement. The K1 mutant model was created by replacing the lysine-76 side-chain in the wild type model with the threonine side-chain using the Richardson library of side chain conformations [43]. The homology model contains errors that are dependent on the percentage of sequence identity between the template and the target and on the number of errors in the template itself. We used several steps to estimate possible errors in our 3D models. The quality of models prepared was assessed using PROCHECK for stereochemical integrity [44], as well as using the Swiss-Model server [45].

Molecular Docking Calculations
Molecular docking calculations were performed using the AUTODOCK package [46,47]. The PfCRT structures (wild type and K76T mutant) were used as receptors. The side chain of residue 76 (LYS for the wild type and THR for the mutant) was kept flexible. The three forms of CQ, viz. neutral, monoprotonated (+1), and bis-protonated (+2), were used as ligands for docking runs. A Lamarckian genetic algorithm (LGA) was used for docking with 100 runs of a population size of 200 for a maximum number of evaluations of 25 million evaluations. The default scoring function set up was used for all calculations. The docked conformations were clustered with an RMSD of 2.5 Å.

Molecular Dynamics (MD) Simulations
All homology models were further refined using molecular dynamics simulations as implemented in Gromacs MD engine [48]. We have used a rhombic dodecahedron box with a minimum distance of 1.0 Å between the protein and box wall. The solvation effect was considered via a TIP3P point charge-based water model. Minimizations were carried out with constraints on all bonds using the P-LINCS algorithm, allowing an integration time step of 2 fs. Long-range electrostatic interactions were calculated using the smooth particle mesh Ewald (PME) method. All short-range nonbonded interactions (including van der Waals terms and the real-space contribution to PME) were truncated at 1.4 nm. A steepest descent minimization algorithm was used, keeping the protein molecule fixed. The energy minimized model was equilibrated in two phases using position restraints applied to all protein heavy atoms (using k pr = 1000 kJ mol -1 nm -2 ). First, a 100 ps canonical equilibration (NVT) was performed using a Berendsen thermostat at 300 K. Subsequently, isothermal-isobaric equilibration (NPT) was performed for 100 ps using an isotropic Berendsen barostat. A production simulation of 10 ns was conducted using an equilibrated system at 300 K without constraints. The amber ff99SB force field was used for all production simulations. Protein-ligand complexes, generated from the top models from molecular docking simulations, were subjected to 1 ns of NVT and NPT equilibration, each in an explicit water medium. These neutral equilibrated systems were subjected to 100 ns production MD simulation without any constraints. The ligands were parameterized with GAFF parameters and AM1-BCC partial atomic charges [49,50]. All trajectory postprocessing and analysis was done using utility tools in the Gromacs package (version 2020.4) and VMD graphics software [51].

Results
Of interest in this work are several features related to structural aspects of the wild type and K1 mutant (Lys76Thr) of PfCRT, which may be of prime importance for CQR. We would like to start discussion of our findings with one important note. The models proposed here do not necessarily address whether PfCRT is a channel or a transporter protein, although the models can serve as starting points toward further assessment of membrane-protein interaction. We have segregated our findings in three major sections. First, the details of the homology models and their quality is presented, followed by molecular dockings of neutral and ionic forms of the CQ on wild-type and mutant protein.
Finally, molecular dynamics simulations of PfCRT-CQ(H/H 2 ) complexes were reported to understand the effect of K76T mutation in chloroquine resistance.

Homology Models of PfCRT
The homology models built using different templates show diverse structural quality (Table 1, Figure 2). In general, multiple template-based models showed the worst stereochemical features, and they were not refined using MD. It should be noted that several of the starting templates had structural clashes, which were transferred in the homology model(s), which prompted us to subject the homology models to MD equilibration runs. The Figure 3 provides top representative homology models. Several key points were to be found amongst various models. The models obtained using template 3RKO and 1EHK had N-and C-terminals at two different ends of the model, whereas the terminals were at the same end for the model obtained using the 1PW4 template. All three models have similar stereochemical features, as shown from Ramachandran plot data. Model 1 is found to have largest helical content amongst the three representative models. Model 1, based on the 1PW4 template, was further used for molecular docking and MD simulations. The choice was made based on the facts that this template belongs to a membrane bound transporter protein and has good overall model quality ( Table 2). model(s), which prompted us to subject the homology models to MD equilibration runs. The Figure 3 provides top representative homology models. Several key points were to be found amongst various models. The models obtained using template 3RKO and 1EHK had N-and C-terminals at two different ends of the model, whereas the terminals were at the same end for the model obtained using the 1PW4 template. All three models have similar stereochemical features, as shown from Ramachandran plot data. Model 1 is found to have largest helical content amongst the three representative models.    Model 1, based on the 1PW4 template, was further used for molecular docking and MD simulations. The choice was made based on the facts that this template belongs to a membrane bound transporter protein and has good overall model quality ( Table 2).

Molecular Docking Calculations
The molecular docking calculations of the CQ (neutral and protonated forms) on the CRT homology models, both wild type and mutant variants, provided support to the CQ-recognition by the mutant version of the PfCRT. The protonated forms of the CQ (+1 and +2) failed to bind to the wild type transporter protein in the site close to the LYS76 residue. These ionic forms are found to be close to the surface of the receptor protein with very unfavorable docking energy (Figure 4). This is due to electrostatic repulsion between the positively charged LYS sidechain and the ionic CQ molecule(s). No such handicap in binding of the ionic forms was observed for the K76T mutant variant. The neutral form of the CQ also showed a significantly larger binding affinity for the mutant variant over the wildtype. The binding energy from the docking calculations is summarized in Table 3. The binding modes are provided in Figure 4. These findings are in line with the increased affinity of the K76T mutant PfCRT to chloroquine in the acidic pH of the DV and explains the role of mutation in CQ resistance.
ing of the ionic forms was observed for the K76T mutant variant. The neutral form of the CQ also showed a significantly larger binding affinity for the mutant variant over the wildtype. The binding energy from the docking calculations is summarized in Table 3. The binding modes are provided in Figure 4. These findings are in line with the increased affinity of the K76T mutant PfCRT to chloroquine in the acidic pH of the DV and explains the role of mutation in CQ resistance.  The 2D diagrams of the binding sites of the CQ drug with the wild type and mutant CRT proteins is provided in Figure 5. The CQ-pharmacophore binds to hydrophobic pockets in both the receptor models. The close contact residues differ between the two proteins, wild type and K76T mutant, although the binding sites are located in the same region on the receptor surface. The CQ binding to the mutant CRT model showed a possible π-stacking interaction involving the quinoline ring and aromatic side chain of phenylalanine-364 residue of the receptor.  The 2D diagrams of the binding sites of the CQ drug with the wild type and mutant CRT proteins is provided in Figure 5. The CQ-pharmacophore binds to hydrophobic pockets in both the receptor models. The close contact residues differ between the two proteins, wild type and K76T mutant, although the binding sites are located in the same region on the receptor surface. The CQ binding to the mutant CRT model showed a possible π-stacking interaction involving the quinoline ring and aromatic side chain of phenylalanine-364 residue of the receptor.

Molecular Dynamics Simulations
For molecular dynamics studies, we have used model PfCRT obtained using glycerol phosphate transporter protein as a template. For all simulations (PfCRT wild and K1 mutant at pH 7.0), the model folds were found to be robust. For all simulations (PfCRT wild and K1 mutant at pH 7.0), the radius of gyration (Rg) and backbone RMSD show little variation, indicating that both models are structurally stable. The amino acid residues are well within the allowed region of Ramachandran plot with few exceptions. A comparison of average protein backbone structure over the 10 ns time frame between K1 and wild type PfCRT model at neutral pH demonstrated no major structural differences (RMSD of 1.3 Å). It is natural as mutating one amino acid (K76T) in PfCRT does not impart significant structural changes; however, the effects of this mutation are rather subtle. Such mutation makes the helix contacting mutation site more rigid, and overall protein structure becomes more compact upon mutation. This reduces electrostatic interaction involving protein and solvent water molecule, which in principle should help drug recognition by providing smaller penalty of desolvation energy of mutant protein upon drug binding, in comparison to the wild type. In order to have qualitative binding affinities of CQ to the mutant and wildtype transporter protein, we have employed MD simulations of the bound CQ complexes. Our calculations showed CQ in a stable bound mode to the K76residue containing the binding site of the mutant protein, whereas for the same binding site in the wild type protein, this binding mode was found to be comparatively weak (Table 4 and Figure 6). The K76T mutation reduces short-range Coulombic interactions between the receptor and CQ ligand and increases the short-range LJ-interactions, providing an overall increase in the affinity of the mutated protein to CQ. The CQ ligand is also found to be more closely associated to the binding site of the mutant protein than the wildtype version.

Molecular Dynamics Simulations
For molecular dynamics studies, we have used model PfCRT obtained using glycerol phosphate transporter protein as a template. For all simulations (PfCRT wild and K1 mutant at pH 7.0), the model folds were found to be robust. For all simulations (PfCRT wild and K1 mutant at pH 7.0), the radius of gyration (Rg) and backbone RMSD show little variation, indicating that both models are structurally stable. The amino acid residues are well within the allowed region of Ramachandran plot with few exceptions. A comparison of average protein backbone structure over the 10 ns time frame between K1 and wild type PfCRT model at neutral pH demonstrated no major structural differences (RMSD of 1.3 Å). It is natural as mutating one amino acid (K76T) in PfCRT does not impart significant structural changes; however, the effects of this mutation are rather subtle. Such mutation makes the helix contacting mutation site more rigid, and overall protein structure becomes more compact upon mutation. This reduces electrostatic interaction involving protein and solvent water molecule, which in principle should help drug recognition by providing smaller penalty of desolvation energy of mutant protein upon drug binding, in comparison to the wild type. In order to have qualitative binding affinities of CQ to the mutant and wildtype transporter protein, we have employed MD simulations of the bound CQ complexes. Our calculations showed CQ in a stable bound mode to the K76-residue containing the binding site of the mutant protein, whereas for the same binding site in the wild type protein, this binding mode was found to be comparatively weak (Table 4 and Figure 6). The K76T mutation reduces short-range Coulombic interactions between the receptor and CQ ligand and increases the short-range LJ-interactions, providing an overall increase in the affinity of the mutated protein to CQ. The CQ ligand is also found to be more closely associated to the binding site of the mutant protein than the wildtype version.  Figure 6. Minimum distance plot between the center of mass of (COM) CQ ligand and the center of masses of LYS76 of wildtype (black) and THR76 of mutant (red) PfCRT from MD production simulation trajectory. The average distance between the COMs of CQ ligand and THR76 residue is smaller than that between that of CQ and LYS76, as observed from MD production simulation.

Discussion
The Lys76Thr mutation in the CRT protein pertains to important molecular recognition to the CQ drug in the Plasmodial parasite. This mutation replaces a lysine side chain that would impart a positively charged local environment with side chain protonation in the acidic pH at DV. Such a local ionic environment hinders CQ binding to the protein.
The CQ molecules will be protonated in the acidic pH of DV. The removal of the lysine by a neutral threonine residue removes the barrier in CQ(protonated)-binding to the CRT protein, which in turn removes the drug from the DV. The molecular docking calculations clearly supports this hypothesis. First, the protonated forms of the CQ drug showed no favorable interactions with the wild type protein model but showed stable binding modes with the K76T mutant model. Second, the molecular dynamics simulations showed higher non-covalent interactions between the neutral CQ-pharmacophore with the mutant protein than the same with the wild type receptors. Together, these two points support the possible enhancement in molecular recognition of the CQ drugs by mutant CRT, thus causing CQ resistance in the malaria parasite.

Conclusions
To summarize, here, we have proposed model structures of wild type and K1-mutant PfCRT protein based on a homology model and molecular dynamics simulations. All the models developed and studied in this report are found to have quite stable folds with well-defined secondary structural elements assigned. The K76T mutation is found to make the helical region more compact around this mutation site. Further, it reduces electrostatic potential of the protein surface. Such changes are more favorable towards recognition of protonated CQ, which should experience more electrostatic repulsion from protonated Lys side chain in wild type protein. This transporter protein is a membrane-bound protein, which separates the highly acidic vacuolar environment from the neutral cytosolic one, and considering different pH effects on the protein structure may provide further key insights involved in CQ resistance. Further studies are in progress to assert the mechanism of chloroquine efflux via PfCRT. Figure 6. Minimum distance plot between the center of mass of (COM) CQ ligand and the center of masses of LYS76 of wildtype (black) and THR76 of mutant (red) PfCRT from MD production simulation trajectory. The average distance between the COMs of CQ ligand and THR76 residue is smaller than that between that of CQ and LYS76, as observed from MD production simulation.

Discussion
The Lys76Thr mutation in the CRT protein pertains to important molecular recognition to the CQ drug in the Plasmodial parasite. This mutation replaces a lysine side chain that would impart a positively charged local environment with side chain protonation in the acidic pH at DV. Such a local ionic environment hinders CQ binding to the protein. The CQ molecules will be protonated in the acidic pH of DV. The removal of the lysine by a neutral threonine residue removes the barrier in CQ(protonated)-binding to the CRT protein, which in turn removes the drug from the DV. The molecular docking calculations clearly supports this hypothesis. First, the protonated forms of the CQ drug showed no favorable interactions with the wild type protein model but showed stable binding modes with the K76T mutant model. Second, the molecular dynamics simulations showed higher non-covalent interactions between the neutral CQ-pharmacophore with the mutant protein than the same with the wild type receptors. Together, these two points support the possible enhancement in molecular recognition of the CQ drugs by mutant CRT, thus causing CQ resistance in the malaria parasite.

Conclusions
To summarize, here, we have proposed model structures of wild type and K1-mutant PfCRT protein based on a homology model and molecular dynamics simulations. All the models developed and studied in this report are found to have quite stable folds with well-defined secondary structural elements assigned. The K76T mutation is found to make the helical region more compact around this mutation site. Further, it reduces electrostatic potential of the protein surface. Such changes are more favorable towards recognition of protonated CQ, which should experience more electrostatic repulsion from protonated Lys side chain in wild type protein. This transporter protein is a membrane-bound protein, which separates the highly acidic vacuolar environment from the neutral cytosolic one, and considering different pH effects on the protein structure may provide further key insights involved in CQ resistance. Further studies are in progress to assert the mechanism of chloroquine efflux via PfCRT.