Insights into the drug resistance induced by the BaDHPS mutations: molecular dynamic simulations and MM/GBSA studies

Dihydropteroate synthase (DHPS) is essential for the folic acid biosynthetic pathway in prokaryotes; the mutation forms for DHPS are found to be relative to the urgent drug resistance problems. In our study, the Bacillus anthracis DHPS (BaDHPS) was selected for molecular dynamics and binding free energy studies to investigate the biochemistry behaviors of the wild-type and mutation form BaDHPS proteins (D184N and K220Q). It is found that the conformational change of the ligand dihydropteroate sulfathiazole binding site in mutation D184N and K220Q systems is mainly attributed from the Loop 1, Loop 2, and Loop 7 regions, and the binding free energy of these mutation systems is lower than that of the wild-type system. Additionally, some important hydrogen bonds of the mutation systems are disrupted during the simulations. But the shortening of the distance between residue Thr67 and the ligand would cause significant change of the binding pose in the K220Q system. These studies of DHPS family will be helpful for further drug resistance investigations. An animated Interactive 3D Complement (I3DC) is available in Proteopedia at http://proteopedia.org/w/Journal:JBSD:24


Introduction
Nowadays, drug resistance has been an urgent problem that severely limits the therapy of current clinical microbial diseases (Boucher et al., 2009;Greenfield & Bronze, 2003). Resistance to sulfonamides is widespread in bacteria and parasites, and generally correlates with mutations to the dihydropteroate synthase (DHPS) gene (Dibbern & Montanaro, 2008;Gibreel & Sköld, 1999;Haasum et al., 2001;Sköld, 2000Sköld, ,2001. DHPS (EC 2.5.1.15) is essential for the folic acid biosynthetic pathway in prokaryotes, in lower eukaryotes such as protozoa and yeast, and in plants (Bermingham & Derrick, 2002;Richey & Brown, 1969). Because humans lack several enzymes for the folate pathway, DHPS has been long exploited as a particularly attractive target for the antimicrobial drugs (Garçon, Bermingham, Lian, & Derrick, 2004;Garçon, Levy, & Derrick, 2006;Goulding et al., 2005;Lawrence et al., 2005;Li et al., 2005). The DHPS family has a classic (α/β) 8 triosephosphate isomerase (TIM) barrel structure, with the active site located at the C-terminal end of the barrel (Achari, 1997;Bermingham & Derrick, 2002). In the DHPS structures, there are two highly conserved loop regions (Loop 1 and Loop 2) which are likely to be of importance to catalysis (Achari, 1997). As a result, many of the drug resistance point mutations are located within the two conserved loops (Baca, Sirawaraporn, Turley, Sirawaraporn, & Hol, 2000;Hampele et al., 1997;Levy, Minnis, & Derrick, 2008;Pemble et al., 2010). Loop 1 and Loop 2 show a high degree of conformational flexibility, causing the inability to observe the DHPS crystal structure determination completely. Then this incompleteness in crystal structure of DHPS will hamper the efforts to understand the structural basis of catalysis mechanism and the interactions between ligand and protein.
Additionally, the Loop 6 and Loop 7 form the wall of binding pocket. These two loop regions are also the major targets for mutation studies. DHPS is encoded by the folP gene, and this enzyme catalyses the pivotal step of folate pathway: the formation of the folate intermediate 7,8-dihydropteroate from p-aminobenzoic acid (pABA) and 6hydroxymethyl-7,8-dihydropterin-pyrophosphate (Richey & Brown, 1969;Shiota, Baugh, Jackson, & Dillard, 1969).
Since 1997, a variety of DHPS crystal structures have been deposited in the Protein Data Bank (PDB) from a few bacterial sources: Staphylococcus aureus, Mycobacterium tuberculosis, Bacillus anthracis, Thermus thermophilus, Escherichia coli, Yersinia pestis, and Streptococcus pneumoniae, and also one from the fungus Saccharomyces cerevisiae (Achari, 1997;Babaoglu, Qi, Lee, & White, 2004;Baca et al., 2000;Hampele et al., 1997;Lawrence et al., 2005;Levy et al., 2008;Yun et al., 2012). B. anthracis served as an important antimicrobial target during the last few years. The crystal structure of B. anthracis dihydropteroate synthase (BaDHPS) was first determined by Babaoglu et al. (2004). Later in 2010 and 2012, some BaDHPS structures complexed with substrates, products, and other inhibitors have been deposited in the PDB (Hevener et al., 2009;Yun et al., 2012). The observation of Achari (1997) and Yun et al. (2012) suggested that some mutations of conserved residues would affect the ligands binding strikingly. Then Babaoglu et al. (2004) found that there was essential water that formed a significant hydrogen bond with the ligand. Recently, it has been observed that the mutations of DHPS will lead to the decrease of drug activity. There is an urgent need to investigate the drug resistance induced by mutations. However, there are no detailed computational studies for the wild-type and mutant forms of BaDHPS. Therefore, to characterize the mechanism of drug resistance and to design new agents according to the mutations are becoming the vital parts of the microbial diseases therapy.
Our current study is aimed to discover the interaction between protein and ligand by using the molecular dynamics (MD) method and molecular mechanics generalized born surface area (MM/GBSA) theory. Moreover, by mutation study, we would explain how the drug resistance mutations affect the affinity of the inhibitor bound to the BaDHPS. As we all know, binding free energy and trajectory analyses have shown to be powerful and valuable tools to understand the mechanisms of the interaction between small molecules and their targets (Gohlke, Kiel, & Case, 2003;Hou, Wang, Li, & Wang, 2011a, 2011b. Here, we report the MD simulation results of BaDHPS, with good agreement with the experimental data, which may be helpful for further drug resistance and de novo drug design investigations.

Initial model of systems
The atomic coordinates and the structural factors for the crystal BaDHPS complex with the inhibitor dihydropteroate sulfathiazole (DHP-STZ) have been deposited in the Protein Data Bank (PDB ID: 3TYE) (Yun et al., 2012), including 273 amino acid per monomer. Chain A (Lys2 to Lys274) of BaDHPS was selected as the initial struc-ture for MD simulations with the ligand DHP-STZ. All the crystal water molecules were retained in the initial model for simulations. Part of the Loop 1 region (Asp31 to Gly36) was invisible in crystal structure 3TYE because of high flexibility, and this part was modeled by the crystal structure of apo BaDHPS (PDB ID: 1TWS) (Babaoglu et al., 2004). Complex systems were generated using PyMOL (DeLano, 2002). In addition, two mutant models of the conserved residues (D184N and K220Q) were prepared by the Swiss PDB viewer software (Guex & Peitsch, 1997). All missing hydrogen atoms of the proteins were added using the LEaP module in the AMBER 11 package (Case et al., 2010). The ff99SB force field was applied to produce the parameters for the protein.
The general Amber force field (Wang, Wolf, Caldwell, Kollman, & Case, 2004) was used to obtain the force field parameters for the inhibitor DHP-STZ, including the Lennard-Jones, torsion, and angle terms. Finally, an appropriate number of sodium counterions were placed in the complex systems, which were solvated in an octahedral periodic box of TIP3P water molecules with a minimum distance of 8 Å between the outermost protein atoms and the walls of the simulation box.

MD simulations
Three systems were prepared for the MD simulations: wild-type BaDHPS/DHP-STZ complex (wt), single-site mutant D184N BaDHPS/DHP-STZ complex (D184N), and single-site mutant K220Q BaDHPS/DHP-STZ complex (K220Q). Two step energy minimizations were performed to all the systems, using SANDER module of AMBER 11 (Pearlman et al., 1995). First, all the water molecules and counterions were minimized with 2000 steps of steepest descent followed by 3000 steps of conjugate gradient. Then the systems were minimized with the same process to remove the bad contacts. The three systems were gradually heated up to 300 K in the canonical (NVT) ensemble, with a velocity of 1 K/ps. Then, they reached equilibration after another 200 ps simulation. Meantime, weak restrains were performed on the Cα atoms of protein and all the atoms of ligands during the first two processes to ensure the accomplishment of the stabilization. Finally, 20 ns MD simulation was carried out with the isothermal-isobaric (NPT) ensemble. SHAKE algorithm (Kräutler, Van Gunsteren, & Hünenberger, 2001) was applied to constrain all bonds involving hydrogen atoms, and particle mesh Ewald method (Darden, York, & Pedersen, 1993) was used to calculate the electrostatic interactions with a cutoff value of 10 Å. The time step was 1 fs. All figures of DHPS in this contribution were generated by Chimera 1.5.3 (Pettersen et al., 2004).

MM/GBSA calculations
The binding free energies of the three enzyme-inhibitor complexes were calculated through the MM/GBSA method (Gohlke et al., 2003;Still, Tempczyk, Hawley, & Hendrickson, 1990). In MM/GBSA, the binding free energy, ΔG bind , is written as the sum of molecular mechanics energy, ΔE MM , salvation free energy contribution to binding, ΔG sol , and the conformational of entropic contribution, ÀTΔS. This can be represented as: where both ΔE MM and ΔG sol can be divided into two parts: The ΔE ele and ΔE vdw in Equation (2) are the electrostatic interaction and van der Waals energy in the gas phase, the ΔG gb and ΔG np in Equation (3) are the polar and nonpolar contributions to the solvation free energy, respectively. The polar contribution of the solvation energy (ΔG gb ) was calculated using the generalized Born model, which parameters were developed by Onufriev, Bashford, and Case (2004) The nonpolar contribution of salvation energy (ΔG np ) was calculated by: The solvent accessible surface area (SASA) in Equation (4) was estimated using ICOSA model in AMBER (Case et al., 2010). And in this study, the values for γ and β were set to 0.0072 kcal mol À1 Å À2 and 0 kcal mol À1 .
The calculations for the binding free energy were accomplished by using the mmpbsa module of AMBER 11 (Case et al., 2010). Then the normal mode analyses were performed to estimate the change of conformational entropy upon the ligand binding (ÀTΔS) via the nmode module of AMBER 11 (Case et al., 2010). For each system, 400 snapshots were extracted from the 20 ns MD trajectory of the complex at an interval of 50 ps. Considering the extremely high computational time, only 100 snapshots of the 400 snapshots were applied for the entropy calculation.
Energy decompositions were performed by the mmpbsa module of AMBER 11 (Case et al., 2010), including the per-residue decomposition and pairwise decomposition. The per-residue decomposition was to separate the energy contribution of each residue from the association of receptor with the ligand into three terms: van der Waals contribution (ΔE vdw ), electrostatic contribution (ΔE ele ), and salvation contribution (ΔG gb + ΔG surf ). And pairwise decomposition calculates the interaction energy between pairs of residues in the system. By analyzing the energy decompositions, we can find out the important interaction within the three systems and calculate the energy.

Cross-correlation analysis
Cross-correlation analysis was used to investigate the extent of correlation motions during the three complex trajectories. Generally, the cross-correlation matrix, C ij , was calculated to reflect the Cα atoms fluctuations relative to their average coordinates. All the 20 ns simulations of the three systems were selected to the calculations by the following equation: In Equation (5), the angle bracket represents an average over the sampled period and Δr i indicates the deviation of the Cα atom of ith residue from its mean position (Amadei, Linssen, & Berendsen, 1993). The value of C ij fluctuates from À1 to 1. Positive C ij values represent a correlated motion between the ith residue and the jth residue, while the negative values describe an anticorrelated motion.
3. Results and discussion 3.1. MD for wild-type and mutation complexes MD simulation of 20 ns was carried out on each system to study the conformational change and to calculate the binding free energy for BaDHPS. All the three systems have reached equilibration after simulation, which prove that the MD simulation method for each complex is reliable during the trajectory. The root mean square deviation (RMSD) values of the three systems with respect to the crystal structure (3TYE) are shown in Figure 1 for the protein and ligand forms, respectively. The backbone RMSD of wt protein is nearly 1.5 Å when compared with the crystal structure. The two mutant proteins, D184N and K220Q, stabilize at around 1.8 Å (Figure 1(a)). In addition, RMSD analyses of the ligand DHP-STZ show similar results. The ligand RMSD values in wt, D184N, and K220Q complexes reach equilibration at 1.0, 1.5, and 2.0 Å (Figure 1(b)). Figure 1(c) shows the comparison between the crystal structure 3TYE and the lowest energy structure of wt form DHPS. It is clear that the MD simulation results of the ligand in wt system are consistent with the experimental results of crystal structure. The location of the essential water in the wt system is the same as that in the crystal structure 3TYE (Yun et al., 2012). In summary, all the mentioned above illustrate that both the mutated protein have a conformational change compared with the wt protein. We should also notice that during the last 5 ns of the trajectories, RMSD values of both protein and ligand in the K220Q system are a bit higher than that in the D184N system.
The important role of the loop regions has been the focus for the investigation of DHPS binding to the ligands. Loop 1 and Loop 2 regions are highly conserved, but Loop 7 shows much lower similarity than Loop 1 and Loop 2 during different DHPS species. Babaoglu et al. (2004) found that one role of Loop 2 was to stabilize the binding of pABA (one of the DHPS substrates) at the active site. After simulations, root mean square fluctuation (RMSF) was measured to investigate the effect of mutation on different areas of the protein during the MD simulations, see Figure 2(a). Three loop regions show large conformational change during the simulations. They are residue Leu26 to Gly37, Gly63 to Val74, and Ile223 to Glu233, which belong to Loop 1, Loop 2, and Loop 7 (Babaoglu et al., 2004), respectively (Figure 2(b)). Residue 184 belongs to the Loop 6 region, but Loop 6 has low RMSF values during the simulations. The three loop regions of the mutation complexes all have larger RMSF values than that of the wt complex. Both the mutations have large RMSF values in the three loop regions; the large conformational change of Loop 1, Loop 2, and Loop 7 in the mutation system leads these regions far away from the binding site (Figure 2(b)). Therefore, we could conjecture that the large conformational change in the three loop regions of mutation complexes would affect the affinity of ligand binding to protein.

Cross-correlation analysis
Cross-correlation analysis was used to study the changes in the internal motions induced by the mutations D184N and K220Q. Overall, it is clear that almost all the interactions of wt protein simulation are filled with positive correlation movements (Figure 3). These movements are shown in red, wine, and yellow colors. By contrast, some internal interactions change into anticorrelation  motions under the influence of mutations. The anticorrelation motions are colored with blue, dark blue, and black, see Figure 3.
In the case of D184N mutation system, it is easy to find out that the anticorrelation motions appear at the C terminal region. Moreover, the residues around the mutation site 184 (most in Loop 6 and α6) in the D184N system have weaker internal interactions than that in the wt system (Figure 3(b)). The crystal structure of BaDHPS (3TYE) reveals that this region is located in the vicinity of ligand DHP-STZ, which is essential to the ligand binding affinity. In the case of K220Q system, the area of anticorrelation motions increases strikingly, not only in the binding site, but also in the other regions (Figure 3 (c)). It is interesting that the mutation K220Q would influence the conformational of the whole protein BaDHPS.
The cross-correlation results suggest that the internal motions of protein BaDHPS are significantly affected by both the mutations, especially the mutation K220Q. The protein tends to be looser in the mutation K220Q system than the D184N system.

MM/GBSA calculation and energy decomposition
To explore the interaction between three kinds of protein BaDHPS (wt, D184N, and K220Q) and the inhibitor DHP-STZ, respectively, the MM/GBSA calculations and energy decompositions were performed using the mmpbsa module of AMBER 11 package. A total of 400 snapshots were extracted from beginning to the end of each system's trajectory for the binding free energy analyses. The predicted binding affinities of the three systems are listed in Table 1.
Overall, as shown in Table 1, the binding free energies of the wt, D184N, and K220Q complexes are À38.77, À34.21, and À24.01 kcal mol À1 , respectively. The wt protein BaDHPS has the highest ability to bind to the ligand DHP-STZ. The binding free energy of mutant D184N decreases a bit, but the K220Q mutant causes the binding ability to decrease dramatically. The binding assay experiments of Yun et al. (2012) and Achari et al. (1997) indicated that both the mutation D184N and the mutation K220Q would decline the affinity of the substrate binding to the BaDHPS. Our computational results also prove the crucial role of the residues Asp184 and Lys220 in the inhibitor binding.
From Table 1 it should be noted that electrostatic interactions play a major role in the binding free energy contributions in the three systems wt, D184N, and K220Q (À103.47, À87.76, and À67.48 kcal mol À1 ), as being compared with the van der Waals interactions (À45.11, À43.53, and À40.12 kcal mol À1 ). In addition, Table 1 also demonstrates that both the two single-site mutants (D184N and K220Q) cause the decline of electrostatic binding free energy, which is the main aspect of binding affinity decrease. In the case of K220Q system, the value of electrostatic binding free energy is almost half as the electrostatic binding free energy of wt system. The nonpolar solvation energies (ΔG np ) are contributed to the ligands binding affinity. The polar solvation energy (ΔG gb ) shows the unfavorable contributions. Consequently, mutations D184N and K220Q would disrupt some electrostatic interactions between protein and ligand.
The residues that contribute to the ligand binding are explored by per-residue binding free energy decomposition analyses. As illustrated in Figure 4, it is obvious that all the residues labeled in the wt system have stronger interactions than the other two mutant systems, especially the residue Asp184 and Lys220. In the case of wt complex, the total binding free energies of these two residues are À8.88 and À4.61 kcal mol À1 , but decreasing to À3.96 and À3.65 kcal mol À1 , À6.10 and À0.03 kcal mol À1 in D184N and K220Q complexes, respectively ( Figure 4). The mutation Asp184 to Asn only influences the interaction between residue 184 and ligand DHP-STZ. In contrast, mutation at the position Lys220 would weaken the interactions or even make  The mutated residue Gln220, colored as magenta in Figure 4 (c), has much lower binding free energy than that residue of the other two systems.
them disappeared, as shown in Figure 4(c). However, the important interaction between Thr67 of K220Q protein and ligand DHP-STZ may induce differences in the binding pose of ligand to protein.
Then the electrostatic energy and van der Waals energy decompositions were analyzed in Figures S1 and S2, which elucidate the residues with binding free energy higher than À1.0 kcal mol À1 . The strong interaction between Lys220 in wt protein and ligand is composed of van der Waals interaction (À2.81 kcal mol À1 ) and electrostatic interaction (À18.67 kcal mol À1 ). In the case of mutation K220Q, both these interaction energies nearly decline to 0. The interaction between Asp184 of wt protein and ligand is totally electrostatic interaction (À13.94 kcal mol À1 ). When the residue Asp is mutated into Asn, weak interaction exists in the D184N complex ( Figure S2). Additionally, it can be observed that the van der Waals interaction between Arg254 and ligand and the electrostatic interaction between Arg234 and ligand decrease in both the mutation systems. The results illustrate that the mutations can disrupt some interactions between the BaDHPS residues and ligand; the residues mainly locate at the loop regions around the ligand binding site (Loop 1, Loop 2, Loop 6, and Loop 7).

Loop regions change under mutations
The hydrophobicity surface of each protein structure combined with ligand DHP-STZ was plotted in Figure 5 in order to find the change of loop regions during each trajectory. Figure 5(a) shows the binding pose of DHP-STZ to the wt protein. The Loop 1, Loop 2, and Loop 7 regions form the wall of hydrophobic tunnel that contains the ligand of protein. This tunnel is long and closed to the ligand for the strong interaction between them. In the case of mutation D184N ( Figure 5(b)), the Loop 7 moves far away from the Loop 1, Loop 2, and the ligand, leading to a looser tunnel than that in the wt system. The interactions between Loop 7 and pentaring, Loop 7 and sulfur atom of sulfathiazole group decrease in the D184N system, which causes the direction of penta-ring changes in this complex. The location of Loop 7 differs dramatically in the K220Q system ( Figure 5(c)). It can be easily found that the volume of binding site increases in K220Q, some of the hydrophobic residues are exposed to the solvent. Additionally, the conserved and essential water molecule for DHPS biochemical reaction is embedded in the binding pocket of both wt and D184N systems. That water approaches the solvent environment directly in the K220Q mutation, which may lead to the decrease of reaction activity.
Moreover, molecular surface module of Chimera software (Pettersen et al., 2004) was used to calculate the SASA of each protein. The SASA calculation results for the three systems are 12,862.45, 12,960.76, and 14,454.80 Å 2 , respectively. Then we compare the SASA values of some residues (Phe189, Pro193, Val208, Leu209, Val226, and Leu234) that contribute to the increase in the mutation systems in Figure S3. The results demonstrate that mutation K220Q will expose some important hydrophobic residues to solvent. As a result, the conserved residue Lys220 is indispensable for keeping the ligand binding pose stable.

Mutations affect the ligand binding pose
The results of the three systems during the MD simulations were analyzed using the ptraj module of the AMBER 11 package. The distance and the angle of hydrogen bonds (Hbonds) calculations were set to 3.5 Å and 120°, respectively. Aiming to find out the reason of the binding pose changes, the H bonds between the inhibitor DHP-STZ and lowest energy structure of wt, D184N and K220Q are represented in Figure 6. Table S1 shows the detail data of these H bonds. Firstly, we focus the analyses on the mutation sites, Asp184 and Lys220. The Asp184 of wt protein can form four strong Hbonds to the ligand DHP-STZ with high occupied percentage and lifetime, but the number of these Hbonds reduces to 1 in the D184N system (Table S1). The Hbonds between the Asp184 of K220Q protein and ligand become weaker than that of wt protein. The residue Lys220 has two Hbonds with the ligand, this number decreases to 1 and 0 in the system D184N and K220Q systems, respectively. Then, under the mutation K220Q, the conformations of the residues around Lys220 change a lot. Figure 6 illustrates that the Hbond between Ser221 and ligand disappears in the K220Q system. Additionally, it should be noted that there is a strong Hbond in the K220Q system between the oxygen atom of Thr67 and the hydrogen atom of the ligand (Figure 6(c)). Then, we measure the distance between the two atoms during the three trajectories, see Figure S4. In the crystal structure 3TYE, this distance is 3.4 Å. This value increases to around 4.0 Å in both wt and D184N systems. The Hbond between the oxygen atom of Thr67 and the hydrogen atom of the ligand exists from beginning to the end of the trajectory of K220Q, with the distance below 2.0 Å. This strong Hbond drags the middle nitrogen atom of ligand to the Loop 2 region as the decrease of distance. In conclusion, the mutation K220Q forms an Hbond with Thr67 and disrupts some Hbonds with protein leading to the change of ligand binding pose.

Conclusion
Drug resistance is the reduction in effectiveness of a drug such as an antimicrobial or an antineoplastic in curing a disease or condition. DHPS is essential for the folic acid biosynthetic pathway in pathogen. Furthermore, some mutation forms of DHPS are relative to the drug resistance.
MD simulations and MM/GBSA calculations were carried out to study the structural changes of BaDHPS and its mutants, that structural changes would be the reason of binding affinity decrease in the mutants D184N and K220Q. Our results demonstrate that large conformational changes of the two mutation systems mainly exist in the Loop 1, Loop 2, and Loop7 regions, which are close to the ligand, and these large conformational changes would cause ligand binding site differing a lot from the wt system. Both the two mutations have lower binding free energies to the ligand, but the binding free energy of K220Q system declines much more than that of D184N. The significant decrease of the binding free energy would pronouncedly influence the binding pose of the ligand. Additionally, some important Hbonds are disrupted in the K220Q system. It should be noted that the Thr67 of K220Q protein forms an Hbond with ligand, which contributes to the conformational change of the binding site and to the binding affinity decrease of the ligand. Finding the cause of ligand affinity decreasing in the mutation proteins will pave the way for further drug resistance investigations on DHPS family.