Protonated Form: The Potent Form of Potassium-Competitive Acid Blockers

Potassium-competitive acid blockers (P-CABs) are highly safe and active drugs targeting H+,K+-ATPase to cure acid-related gastric diseases. In this study, we for the first time investigate the interaction mechanism between the protonated form of P-CABs and human H+,K+-ATPase using homology modeling, molecular docking, molecular dynamics and binding free energy calculation methods. The results explain why P-CABs have higher activities with higher pKa values or at lower pH. With positive charge, the protonated forms of P-CABs have more competitive advantage to block potassium ion into luminal channel and to bind with H+,K+-ATPase via electrostatic interactions. The binding affinity of the protonated form is more favorable than that of the neutral P-CABs. In particular, Asp139 should be a very important binding site for the protonated form of P-CABs through hydrogen bonds and electrostatic interactions. These findings could promote the rational design of novel P-CABs.


Introduction
The gastric H + ,K + -ATPase (proton pump), primarily responsible for gastric acid secretion, is the key therapeutic target for the ulcer diseases such as duodenal ulcers, gastric ulcers, gastro esophageal reflux disease (GERD), Zollinger-Ellison syndrome (Z-E), and gastritis [1][2][3]. As a member of the P 2 -type ATPase family, H + ,K + -ATPase is a dimeric heterodimer composed of a subunit of about 1035 amino acids with 10 transmembrane (TM) segments and b-subunit glycoprotein with 290 amino acids [4,5]. By cyclic phosphorylation and dephosphorylation of the catalytic subunit, H + ,K + -ATPase undergoes conformational changes between E 1 and E 2 . With phosphorylation, H + ,K + -ATPase E 1 conformation binding hydronium ion E 1 PNH 3 O + changes to E 2 PNH 3 O + form. After release of H 3 O + and binding of K + on the extracytoplasmic surface, the E 2 PNK + conformation is formed and then converts to the E 1 K conformation with the dephosphorylation. The E 1 K conformation releases K + to the cytoplasmic side, then rebinding of H 3 O + occurs to complete the transition cycle [6]. The H + , K + -ATPase engages in 2K + /2H + /1ATP electroneutral ion exchange, generating a million-fold H + -gradient across the mammalian canalicular membrane of the parietal cell [7,8].
So far the structural investigations of H + ,K + -ATPase have lagged behind the pharmacological studies. The structure of human gastric H + ,K + -ATPase is unknown, and the structure of pig gastric H + ,K + -ATPase is poorly defined, being currently limited to a resolution of 7 Å (PDB code: 3IXZ [32], resolution: 6.5 Å ; PDB code: 2XZB [33], resolution: 7 Å ). So the aim of this study is to  model H + ,K + -ATPase structure by homology modeling and for the first time to investigate the interactions between the protonated form of P-CABs and H + ,K + -ATPase using molecular docking, molecular dynamics and MM/GBSA calculation methods.

Homology modeling
The sequences of the human (1035 amino acids) and pig (1033 amino acids) gastric H + ,K + -ATPase receptor were taken from the Swiss-Prot Database (ID: P20648 and P09626) [34]. Through BLAST online method (http:// blast.ncbi.nlm.nih.gov/Blast.cgi) [35] from the Protein Data Bank [36], the crystal structure of Na + ,K + -ATPase in the E 2 P state (PDB code: 2ZXE) [37] was used as a template. The sequence alignment was performed with the ClustalW2 algorithm [38]. The homology model of human H + ,K + -ATPase was generated using MODELLER9v4 [39]. The resultant structure of the H + ,K + -ATPase was subject to the Protein Preparation Wizard module in Schrödinger [40] as follows: adding hydrogens, assigning partial charges, and minimizing using the OPLS-2005 force field [41] until RMSD 0.30 Å . The final optimized model was validated using the program PROCHECK [42] to assess the quality of the stereochemistry of the protein structure.

Preparation of ligands
LigPrep of Schrödinger software suit [43] was used for the preparation of ligands: generating 3D structures from 2D (SDF) representation, and performing a geometry minimization. The ligands were subjected to energy minimization using MacroModel module of Schrödinger with Merck Molecular Force Field (MMFFs). Truncated Newton Conjugate Gradient (TNCG) minimization method was used with 500 iterations and convergence threshold of 0.05 (kJ/mol). While Epik [44] was used to generate possible ionization states at pH 7.061.0.

Molecular docking
The docking simulations were performed using the software Glide (XP mode) [45,46]. Previous biochemical and mutagenesis studies [18,33,[47][48][49][50] suggest that Cys813 in pig H + ,K + -ATPase (corresponding to Cys815 in human H + ,K + -ATPase) is the key amino acid residue in the luminal cavity. Therefore, dimensions for the cubic boundary box centered on the centroid of Cys815 were set to 20 Å620 Å620 Å . The scaling factor for protein van der Waals radii was 1.0 in the receptor grid generation. After docking calculations, for each ligand, the best pose was chosen and scored using the proprietary GlideScore function.

QM/MM optimization
The final docking complexes were energetically optimized by QM/MM method. QM/MM calculations were carried out using the QSite program [51,52] of the Schrödinger suite. The ligands were defined as QM region calculated by the density functional theory DFT/B3LYP (6-31G* basis set). The receptor as MM region was minimized with Truncated Newton algorithm (maximum cycles as 1000; gradient criterion as 0.01). The OPLS 2005 all-atom force field was employed.

Molecular dynamics
The optimized docking models were subjected to molecular dynamics simulations using Desmond [53,54]. The system was embedded in a fully hydrated POPC (1-palmitoyl-2oleoyl-snglycero-3-phosphatidylchlorine) bilayer membrane and solvated with an orthorhombic box of SPC water molecules (buffer distance: 10 Å610 Å612 Å ). Counter-ions (Na + ) were added to neutralize the system and a physiological salt concentration of 0.15 M NaCl was introduced. The final system was composed of approximately 130,000 atoms. Each model was equilibrated in MD for 10 ns. Then two K + were added in the luminal domain near channel (close to Asp134) replacing two Na + , and potassium ion competition molecular dynamics simulations were carried out for 10 ns. Before the simulation, the models were relaxed as follows: (1) two minimization steps (restraining the solute and unrestrained minimization) with maximum runs of 2000 and the convergence threshold for minimization set to 1 kcal/mol/Å . The minimization method was a hybrid of the steepest decent and limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithms; (2) after minimization, the simulation in the NVT ensemble was run restraining all solute heavy atoms with temperature of 10 K for 20 ps, using Berendsen thermostat; (3) a simulation in the NPT ensemble restraining all solute heavy atoms with temperature of 10 K and 300K for 50 ps, respectively; (4) a simulation in the NPT ensemble, no restraints, with temperature of 300 K and simulation time of 200 ps. Then 10 ns MD production runs (time step: 2.0 fs) were performed through NPT ensemble at 300 K with 1.0132 bar pressure.
Smooth particle mesh Ewald method (Ewald tolerance: 1e-09) was employed to treat the long-range electrostatic interactions and a 9 Å radius cut off was used for coulombic short range interactions. The energies and frames of each trajectory were recorded every 1 ps and 5 ps, respectively. After 10 ns MD simulations with two K + , SCH28080 and TAK438-1 systems with H + ,K + -ATPase and POPC membrane in 3M KCl resolution were equilibrated in MD for 40 ns and then run 100 ns disassociation molecular dynamics, respectively. MD trajectory analysis was performed using Desmond utilities and VMD [55]. The ligand-protein complexes were visualized using PyMOL [56] and analyzed with Ligand Interactions module embedded in Maestro 9.3 [57].

MM/GBSA calculations
For each system, binding free energy (DG bind ) calculations were performed for 100 snapshots extracted from the last 1 ns stable MD trajectory using molecular mechanics-generalized Born surface area (MM/GBSA) method. MM/GBSA procedure in Prime program [58,59] was used to calculate DG bind of the docked inhibitors according to the following equations [60]: Where DE MM is the difference of the gas phase MM energy between the complex and the sum of the energies of the protein and inhibitor, and includes DE internal (bond, angle, and dihedral energies), DE Elect (electrostatic), and DE VDW (van der Waals) energies. DG solv is the change of the solvation free energy upon binding, and includes the electrostatic solvation free energy DG GB (polar contribution calculated using generalized Born model), and the nonelectrostatic solvation component DG SA (nonpolar contribution estimated by solvent accessible surface area). TDS is the change of the conformational entropy upon binding. This term was calculated using normal-mode analysis Rigid Rotor Harmonic Oscillator (RRHO) contained in MacroModel module [61]. DG ' bind neglects the effect of entropy contributions, while DG bind includes contributions from loss of ligand translational, rotational and vibrational entropy (TDS).

Results and Discussion
The three-dimensional structure of Na + ,K + -ATPase in the E 2 P state (PDB code: 2ZXE; resolution: 2.4 Å ) [37] was selected as a template, which share 63.9% identity and 78.7% similarity to human H + ,K + -ATPase on the basis of sequence alignment analysis ( Figure 2). The human gastric H + ,K + -ATPase model is shown in Figure 3. It comprises ten transmembrane helices (TM1-TM10), in which binding sites are located, and three cytoplasmic domains: the nucleotide-binding (N), phosphorylation (P) and actuator (A) domains [33]. The stereochemistry of the homology model was assessed using Ramachandran plot, which indicates that 93.8% of the residues were located in the most favored zones, 6.1% in allowed regions, 0.0% in generously allowed regions and 0.1% in disallowed regions ( Figure S1). The dihedrals, covalent and overall G-factors of this model are -0.17, 0.01 and -0.08, respectively. The PROCHECK G-factor should be above -0.5 ideally for the homology model to have a good acceptance score.

Molecular dynamics simulations
After molecular docking and QM/MM optimization (Glide docking Gscore, Glide energy and QM/MM energy were listed in Table S1), MD simulations for ten P-CAB-H + ,K + -ATPase complexes with two K + in NaCl aqueous solution were run for a duration of 10 ns. To explore the dynamic stability of complexes and to ensure the rationality of the sampling method, root-meansquare deviations (RMSD) for the backbone atoms from the starting structure were analyzed, as shown in Figure 4. As can be seen in the plots, after 8 ns, the RMSD of each system tends to converge, indicating that the systems are stable and equilibrated. During the whole MD simulation process, the potassium ions were blocked into the luminal channel by all P-CABs, which acted just like goalkeepers. The distances between one of potassium ions (K1) and the nitrogen atom (N15475) of -CN in SCH28080 (or nitrogen atom (N15494) of -NH 2 + CH 3 in TAK-438-1) are shown in Figure 5 with the interaction modes at 0 ns and the time of the shortest distances (SCH28080: 7.570 ns, distance 2.744 Å ; TAK-438-1: 3.180 ns, distance 6.218 Å ). Due to the electrostatic repulsive interaction, the protonated form (+1 charged) can block K + more easily than the neutral P-CABs. From 4.218 ns, the  distance between K + and TAK-438-1 is larger than 20 Å , while this time point is 9.421 ns for SCH28080 system ( Figure 5A). Furthermore, analyses of root-mean-square fluctuation (RMSF) versus the residue number for SCH28080 and TAK-438-1complexes are illustrated in Figure 6. All complexes share similar RMSF distributions and similar trends of dynamic features. They are relatively rigid in the active site region (residues Leu143 in TM2, Ala337 in TM4, Tyr801 in TM5, Leu811 in the TM5-6 loop, and Cys815 in TM6) as reported in the literatures 18,33,[47][48][49][50] (residue IDs in pig H + ,K + -ATPase correspondingly are Leu141, Ala335, Tyr799, Leu809, and Cys813). The fluctuation of the active site for TAK-438-1 complex (except for Leu143) is smaller than SCH28080 complex. RMSF values of the important amino acids in different P-CAB system are listed in Table 2.

Interaction modes of P-CABs with H + ,K + -ATPase
To investigate P-CAB interactions in the binding site, the average structures from last 1 ns MD trajectory were compared ( Figure 7). After MD simulation, the pose of SCH28080 was very similar to that reported by Abe et al. [33]. The binding site is in the luminal cavity, surrounded by Ala337, Tyr801 and Cys815, which have hydrophobic interactions with SCH28080 ( Figure 7A). The hydrogen bond (H-bond) interactions, which play an important role in P-CABs binding to H + ,K + -ATPase, are listed in Table 3 and shown in Figure 7. Oxygen atom of SCH28080 and nitrogen atom of AZD0865 formed hydrogen bonds with Tyr930. TAK-438 was found to form H-bond with Tyr801, and has p-p stacking interaction with Phe334 and p-cation interaction with Arg330. The protonated forms have different interaction modes compared to the neutral P-CABs. There are strong H-bond interactions between Asp139 and TAK-438-1 (distance 1.72760.012 Å ) and Soraprazan-1 (distance 1.62360.006 Å and 2.13360.015 Å ). Asp139 also has negative charged interactions with all protonated P-CABs (+1 charged). AZD0865-1 interacts with Asp138 via two H-bonds (distance 1.69560.009 Å and 1.90360.013 Å ). Both Soraprazan-1 and Revaprazan-1 have H-bonds with Phe990. TAK-438-1 formed hydrogen bonds with Asp139, Tyr927 and Asn991 simultaneously.

Binding free energy of P-CABs
The binding free energies of all ten systems were calculated by MM/GBSA method. As listed in Table 4, the DG ' bind and DG bind values all show that the binding affinity of protonated P-CABs is more favorable than that of neutral P-CABs. The order of favorable binding interaction is TAK-438-1.Soraprazan-1.Revaprazan-1. SCH28080-1.AZD0865-1.TAK-438.AZD0865.SCH28080< Soraprazan.Revaprazan, which is consistent with the experimental results at different pH. According to the protonated form percentage of P-CABs in Table 1, the binding free energies (DG bind ) of protonated and neutral mixtures of P-CABs at different pH values were calculated. The correlation coefficient R between DG bind of the mixture and pIC 50 (negative logarithms of 50% inhibition concentration) is 20.837 (Figure 8).
To estimate which energy term has most impact on the binding affinities, the four individual energy components (DE Elect , DE VDW , DG GB , and DG SA ) were carefully compared. From Table 4, it can be seen that both the van der Waals (DE VDW ) and the electrostatic (DE Elect ) contributions are essential for P-CABs binding to H + , K + -ATPase. For the neutral form, the contributions of DE VDW are more favorable than DE Elect term. But for protonated P-CABs, the major favorable contributor is DE Elect term, which is far larger than for the neutral form. Although polar solvation (DG GB ) term of protonated form is also far larger than for the neutral form, which could offset DE Elect , the net electrostatic contributions (DE Elect + DG GB ) are more favorable for the protonated compounds. Thus, the electrostatic contribution is more crucial for distinguishing the binding affinities among P-CAB complexes. Except for Soraprazan-1 and AZD0865-1, nonpolar solvation term (DG SA ) contribute favorably to binding. Among all P-CABs, DE Elect and DG SA of TAK-438-1 is the most favorable (2145.8861.80 and 29.6160.30 kcal/mol). In addition, the contributions of the conformational entropy (TDS) are very similar for the protonated compounds (22.3460.21,23.3260.29 kcal/mol).
The binding free energy between the protein and P-CABs was decomposed into the contribution of each residue, which provides quantitative information of the key residues related to the detailed binding mechanism. From the energy comparison of residues in binding sites (Figure 9), it can be observed that there is the distinct difference between protonated and neutral forms of P-CABs. The energies of residues Tyr142, Tyr804, Leu811, Leu813, Gly814, and Tyr930 interacting with SCH28080-1 (20.4060 (Figure 9A), which in particular include the key residues Leu143, Tyr801, and Cys815 (all above 25.00 kcal/mol for SCH28080-1). This might be the reason that the activity of SCH28080 at pH 6.5 is dramatically higher than that at pH 7.5. Due to strong hydrogen bond and electrostatic interactions, the energy contributions of Asp139 to protonated form are all more favorable than those to neutral P-CABs (214.6360.51 kcal/ mol to Soraprazan-1, -11.2960.35 kcal/mol to Revaprazan-1, 29.5760.26 kcal/mol to AZD0865-1). In particular, the highest value (223.2860.64 kcal/mol) is reached for TAK-438-1. Even though less favorable than their protonated forms, the interactions with Asp139 are also strong to Soraprazan (25.6960.77 kcal/mol) by H-bond and to AZD0865 (27.1460.41 kcal/mol) by electrostatic interactions. Asp139 could be the key residue in the binding sites of P-CABs, which was not paid attention to previously. TAK-438 has more favorable interactions with the key residues Val333, Tyr801, Leu811, and Cys815 (24.8860.13, 23.5960. 16, 22.8660.09, 22.8460.17 kcal/mol). However, the binding free energy of TAK-438-1 is higher than TAK-438 because TAK-438-1 has strong

Disassociation molecular dynamics
After 100 ns molecular dynamics (RMSD for the backbone atoms and RMSF for the residues of SCH28080 and TAK-438-1complexes are illustrated in Figure S2 and S3), SCH28080 occurred partial disassociation from H + ,K + -ATPase and DG bind of SCH28080 decreased to 227.8860.74 kcal/mol (Table 4). Because water molecules have went into the ion channel, the interactions between SCH28080 and residues Val333 to Val340 in TM4 were cut off ( Figure 10A and Figure 11). In contrast to SCH28080, DG bind of TAK-438-1 after 100 ns MD was 266.486 1.87 kcal/mol (Table 4) and similar to that in 10 ns MD. TAK-438-1 still interacted with Asp139 (225.4860.90 kcal/mol), Tyr927 (211.8161.28 kcal/mol) and Asn991 (29.336 0.31 kcal/mol) by H-bonds and blocked water molecules into the ion channel ( Figure 10B and Figure 11). Therefore TAK-438-1 could have the long dwell time and disassociate from H + , K + -ATPase very slowly, which make it a likely competitor for PPIs.

Conclusions
In the present study, to clarify the interaction mechanism between protonated P-CABs and H + ,K + -ATPase, MD simulations and MM/GBSA binding free energy calculations were conducted. There is a distinct difference in interaction mode between protonated and neutral P-CABs. With positive charges on nitrogen atoms, protonated forms of P-CABs have a competitive advantage in blocking potassium ion into the luminal channel and binding to H + ,K + -ATPase mostly via electrostatic interactions. According to the binding free energy analysis, the binding affinity of protonated Table 4. The binding free energies of P-CAB complexes (kcal/mol). forms is more favorable than that of neutral P-CABs. Asp139 in particular should be a very important binding site for protonated forms of P-CABs through hydrogen bonds and electrostatic interactions. The energy contributions of Asp139 for protonated forms are all more favorable than those for neutral P-CABs. Thus, protonated form is the potent form of P-CABs. To enhance P-CABs pKa values and binding with the key residue Asp139 would help increase the inhibition activity. The 100 ns disassociation molecular dynamics suggest that TAK-438 in protonated form could have the long dwell time and disassociate from H + , K + -ATPase very slowly, which make it a likely competitor for PPIs. The findings in this work provide a better structural understanding of the binding sites in H + ,K + -ATPase and a basis for further rational design of novel P-CABs. Figure S1 Ramachandran plot of the human H + ,K + -ATPase model. (TIF) Figure S2 RMSD for the backbone atoms of the SCH28080 and TAK-438-1complexes in 100 ns disassociation molecular dynamics.