- 1National Supercomputing Center in Zhengzhou, Zhengzhou University, Zhengzhou, China
- 2School of International Studies, Zhengzhou University, Zhengzhou, China
- 3Department of General Surgery, Affiliated Cancer Hospitalof Zhengzhou University, Henan Cancer Hospital, Zhengzhou, China
Gastrin plays important role in stimulating the initiation and development of many gastrointestinal diseases through interacting with the cholecystokinin 2 receptor (CCK2R). The smallest bioactive unit of gastrin activating CCK2R is the C-terminal tetrapeptide capped with an indispensable amide end. Understanding the mechanism of this smallest bioactive unit interacting with CCK2R on a molecular basis could provide significant insights for designing CCK2R antagonists, which can be used to treat gastrin-related diseases. To this end, we performed extensive Gaussian accelerated molecular dynamics simulations to investigate the interaction between gastrin C-terminal pentapeptide capped with/without amide end and CCK2R. The amide cap influences the binding modes of the pentapeptide with CCK2R by weakening the electrostatic attractions between the C-terminus of the pentapeptide and basic residues near the extracellular domain in CCK2R. The C-terminus with the amide cap penetrates into the transmembrane domain of CCK2R while floating at the extracellular domain without the amide cap. Different binding modes induced different conformational dynamics of CCK2R. Residue pairs in CCK2R had stronger correlated motions when binding with the amidated pentapeptide. Key residues and interactions important for CCK2R binding with the amidated pentagastrin were also identified. Our results provide molecular insights into the determinants of the bioactive unit of gastrin activating CCK2R, which would be of great help for the design of CCK2R antagonists.
1 Introduction
Gastrin is a peptide hormone that can regulate gastric acid secretion and stimulate gastric mucosal growth (Dockray, 2004; Dimaline and Varro, 2014). Many variants of gastrin with various amino acid lengths are existing in the human body and the main forms contain 17 or 34 amino acid residues with an amidated C-terminus (amidated gastrins) (Calatayud et al., 2010; Dimaline and Varro, 2014). Amidated gastrins with different amino acid lengths have full biological activities but varied circulating half-lives, while non-amidated forms have almost no biological activities (Calatayud et al., 2010; Dimaline and Varro, 2014). The essential bioactive unit of gastrin is confirmed as the C-terminal tetrapeptide amide sequence (Trp-Met-Asp-Phe-NH2) (Dockray, 2004). Gastrin exerts its biological functions via binding to cholecystokinin 2 receptor (CCK2R), a member of class A G protein-coupled receptor (GPCR) superfamily. CCK2R has another endogenous ligand, CCK, which shares the same C-terminal five residues and amide end as gastrin and is responsible for pancreatic enzyme secretion (Calatayud et al., 2010; Dimaline and Varro, 2014). CCK and gastrin are the only family members of regulatory peptides in mammals (Rehfeld, 2017). The bioactive unit of CCK is its C-terminal seven amino acids amide sequence (Rehfeld, 2017). Apart from the normal physiology effects of gastrin activating the CCK2R signaling pathway, it is also found that CCK2R can be activated by gastrin in many gastrointestinal diseases including peptic ulcer disease, malignant tumors of the colon, pancreatic and gastric cancers, resulting in cell proliferation (Fourmy et al., 2011; Dimaline and Varro, 2014; Maddalo et al., 2014; Smith et al., 2016). Among these diseases, gastric and pancreatic cancers are the third and seventh leading causes of cancer deaths, respectively (Bray et al., 2018). Interrupting the interaction between gastrin and CCK2R should be feasible to treat these deadly cancers. However, despite many strategies that have been explored to develop CCK2R antagonists, there are still no approved drugs. At present, surgical resection is the only effective therapy to treat these malignancies (Gilliam and Watson, 2007; Hussain, 2016; Bray et al., 2018; Dolcetti et al., 2018; Sitarz et al., 2018). Understanding the mechanism of gastrin interacting with CCK2R as well as the dominant role of the C-terminal amide cap should give significant insights into developing CCK2R antagonists.
Considerable experiments have been performed to explore the activation mechanism of CCK2R in the past 30 decades (Silvente-Poirot and Wank, 1996; Bläker et al., 1998; Silvente-Poirot et al., 1999; Pannequin et al., 2002; Irina et al., 2007; Stone et al., 2007; Foucaud et al., 2008; Song et al., 2013; Ritler et al., 2019). Methods include modifying residues in CCK2R, gastrin and CCK to identify the key residues accounting for interaction. It is found that gastrin with sulfated or non-sulfated tyrosine has similar binding affinities with CCK2R (Dockray, 2004; Dimaline and Varro, 2014). The negatively charged pentaglutamic acid sequence in the gastrin is important for gastrin’s stability in human blood plasma and can increase its binding affinity with CCK2R. Segment replacement combined with site-directed mutagenesis identified five residues (Q204CVHRW209) in the second extracellular loop of rat CCK2R playing major roles in the selectivity of gastrin (Silvente-Poirot and Wank, 1996). Point mutations also found some residues in CCK2R can significantly alter the activity of gastrin and CCK(Silvente-Poirot et al., 1998; Anders et al., 1999; Bläker et al., 2000; Galés et al., 2003; Irina et al., 2007; Foucaud et al., 2008; Willard et al., 2012; Magnan et al., 2013; Lipiński et al., 2018). Very recently, three cryo-EM structures of CCK2R in complex with gastrin-17 and Gq or Gi2 (PDB codes: 7F8W, 7XOW and 7F8V) were obtained, which presented the complex conformations at atomic level (Zhang et al., 2021; Ding et al., 2022). Besides experiment methods performed to identify key sites contributing to gastrin/CCK-CCK2R interactions, computation methods including docking and molecular dynamics (MD) simulations were also performed to delineate the complete binding pose and analyze the main interactions (Langer et al., 2005; Marco et al., 2007; Stone et al., 2007; Foucaud et al., 2008; Magnan et al., 2013). However, these MD results are limited by the simulation time of tens of ns which may not be sufficient for conformational adjustment from the initial docking pose. Besides, the molecular dissection of the indispensable role of the amide cap still remains elusive.
Benefiting from the encouraging development of computational power and computational methods (Miao and McCammon, 2016; Mori et al., 2016; Yang et al., 2016; Eastman et al., 2017; Miao and McCammon, 2017; Yang et al., 2019), the whole processes of protein folding, biomolecular large-scale conformational transitions and biomolecular recognition in the explicit solvent can be simulated within an acceptable time (Lindorff-Larsen et al., 2011; Eastman et al., 2017; Miao and McCammon, 2017; Palermo et al., 2017). In the present work, we applied Gaussian accelerated MD (GaMD) simulations (Pang et al., 2017; Bhattarai et al., 2020), an advanced conformational sampling method, to explore the influence of the amide cap on the conformational dynamics of gastrin C-terminal pentapeptide as well as its binding modes with CCK2R. We showed the molecular basis of the gastrin bioactive unit interacting with CCK2R and the mechanism of the indispensable role of the amide cap in gastrin fulfilling its physiological functions. These results should provide deep understanding of the activation mechanism of CCK2R and practical guides for designing CCK2R antagonists.
2 Materials and methods
2.1 System preparation
CHARMM-GUI (Jo et al., 2008) was used to prepare the simulation systems from the cryo-EM structure of CCK2R in complex with gastrin-17 (PDB code: 7F8W) (Zhang et al., 2021). First, the bioactive unit structure of C-terminal pentagastrin (Gly-Trp-Met-Asp-Phe) was truncated from the structure of gastrin-17 and the amidated (termed as G5NH2) and non-amidated (termed as G5) C-termini were built. The neutral N-terminus was used for G5/G5NH2 to imitate the primitive form of pentagastrin in gastrin-17 and minimize the influence of the N-terminus on the binding modes. Besides, the disulfide bond between Cys1273.25 and Cys205ECL2 (superscript indicates nomenclature according to Ballesteros—Weinstein numbering system (Isberg et al., 2015)) in CCK2R was also built. Then the Membrane Builder plugin (Wu et al., 2014) was used to pack the explicit POPC lipid bilayer around CCK2R-G5 and CCK2R-G5NH2 complexes, respectively. Subsequently, 20 Å thickness TIP3P water molecules were placed above and below the lipid bilayer. The protonation states of all residues in G5, G5NH2 and CCK2R were set at PH 7. Additional Na+ and Cl− ions were added to all simulation systems to keep the environment neutral and the salt concentration at 0.15 M.
2.2 Simulation protocols
All simulations were performed by NAMD v2.13 (Phillips et al., 2005) with GPU acceleration. Protein and lipid atoms were parameterized with CHARMM36m force field (Huang et al., 2017). Periodic boundary conditions were applied in all three directions. Short-range non-bonded interactions were computed every step using a cutoff of 12 Å with a switch distance of 10 Å. Long-range electrostatic interactions were calculated every two steps using the Particle-Mesh Ewald (PME) (Essmann et al., 1995) method with a grid spacing of 1 Å−1. The integration time step was set to 2 fs with all bonds involving hydrogen atoms constrained by the SHAKE algorithm (Weinbach and Elber, 2005). Constant temperature of 300 K and constant pressure of 1 atm were kept by langevin dynamics and Nosé-Hoover Langevin piston method (Feller et al., 1995), respectively.
All systems were firstly undergone 5,000 steps of conjugate gradient energy minimization to remove the steric clashes in the system. Then 300 ps gradually heating process to rise the temperature from 0 to 300 K in the NVT-ensemble and 1 ns NPT-ensemble MD simulations used to adjust the volume of each system were performed. In the last two steps, harmonic potential with the force constant of 1.0 kcal/mol/Å2 was applied to all non-hydrogen atoms of solutes to minimize their structures change. Another 1 ns NPT-ensemble MD simulations were performed with the harmonic potential gradually decreased to 0 kcal/mol. Finally, GaMD simulations were carried out for thoroughly conformational sampling.
GaMD is an unconstrained enhanced conformational sampling technique which works by adding a harmonic boost potential to reduce the energy barriers at the system’s potential surface. When the system potential energy (V) is lower than the threshold energy E, a boost potential (ΔV) is added and the system potential is modified as.
where
2.3 Trajectory analysis
Based on the GaMD production simulations, the 2D potential of mean force (PMF) profiles were calculated using PyReweighting toolkit (Miao et al., 2014) with the reweighting method of cumulant expansion to the second order. The two parallel GaMD trajectories of each system were combined first before reweighting. Two reaction coordinates were defined as distances between Cα atoms in N/C-terminal residues in G5/G5NH2 and the geometrical center of the transmembrane domain (calculated for backbone atoms) of CCK2R. The residues ID of each transmembrane domain in CCK2R is listed in Supplementary Table S1. Principal component analysis (PCA) was performed using the Wordom software (Seeber et al., 2011) to differentiate the dominant conformational features of CCK2R induced by G5 and G5NH2 binding. The correlations of the atomic fluctuations in the CCK2R-G5/G5NH2 complex are calculated by the Linear Mutual Information (LMI) (Lange and Grubmüller, 2006) method within the Bio3D package (Grant et al., 2021).
Residues contributing to CCK2R-G5/G5NH2 binding were denoted by MM-PBSA method with MMPBSA.py program (Miller et al., 2012):
where
The binding sites are defined as the atoms of CCK2R located within 4.0 Å of the heavy atoms in G5/G5NH2. In the MD trajectory frames, the occurrence rates of each heavy atom located in the binding sites were calculated, and one residue is denoted as in the binding sites only if the sum of the occurrence rates of every heavy atom in it is larger than 1. Key interactions of hydrogen bonds and salt bridges (Yang et al., 2018) formed between CCK2R and G5/G5NH2 are also analyzed. The hydrogen bond is detected by VMD according to the criteria of the distance of donor-acceptor less than 3.5 Å and the angle of donor-H-acceptor larger than 150°. The salt bridge is defined as the distance between the oxygen atom of the carboxyl group and the nitrogen atom at the side chain of basic residues (Arg and Lys) is less than 4.0 Å. Besides, the occurrence rates of hydrogen bonds and salt bridges in the trajectory frames should be larger than 0.5. Residues which locate in the binding site and made binding free energy contributions larger than 1 kcal/mol are denoted as key residues for binding.
3 Results and discussion
3.1 Free energy profiles of CCK2R-G5 and CCK2R-G5NH2
To clarify the binding mechanism of the smallest bioactive unit of gastrin to CCK2R and why the C-terminal amide cap is indispensable for gastrin activating CCK2R, we performed a total of 2 μs (two individual 1 μs trajectories) GaMD production simulations on the CCK2R-G5 and CCK2R-G5NH2 complexes, respectively. From the GaMD production simulations, 2D PMF profiles along the reaction coordinates of distances between the N/C-terminus of G5/G5NH2 (Cα atoms of G13 and F17, respectively) and the geometrical center of the TM domain of CCK2R were calculated and presented in Supplementary Figures S1, S2 and Figure 1. Supplementary Figure S1 shows that the lowest free energy regions with PMF less than 1 kcal/mol are very similar from 800 ns to 1 μs (combined two parallel trajectories), which suggests the convergence of present GaMD simulations. Besides, the lowest PMF regions for CCK2R-G5NH2 complexes are around the initial site in the cryo-EM conformation, while having some distance from the initial site for CCK2R-G5 complexes (Figures 1A, B). Supplementary Figure S2 shows the PMF profiles calculated from the individual trajectory. From Supplementary Figures S2A, B, we can find that different minimal PMF regions were obtained from different trajectories for the CCK2R-G5 system. To find why the PMF files calculated from individual trajectories differed, we depicted the time evolutions of two reaction coordinates in each trajectory in Supplementary Figure S3. From Supplementary Figure S3A, we can find that for the reaction coordinate of TM center—G13 distance, a significant difference occurred after 150 ns between two trajectories, which indicates different conformational spaces sampled by two trajectories. However, both reaction coordinates have stable fluctuations after 400 ns (Supplementary Figures S3A, B). PMF profiles calculated from merged trajectories (Figure 1A) showed that the minimal PMF region obtained in trajectory 1 (Supplementary Figure S2A) was lower than that obtained in trajectory 2 (Supplementary Figure S2B). To further verify the minimal PMF region found by trajectory 1, we performed another 300 ns simulations (termed trajectory 3) from the conformational state at 700 ns of trajectory 1 with different randomized initial atomic velocities. The calculated PMF profile is depicted in Supplementary Figure S2C, and the time evolutions of two reaction coordinates are depicted in Supplementary Figure S3. We can find that the minimal PMF region and the time evolutions of two reaction coordinates in trajectory 3 are much similar to those in trajectory 1. Adding trajectory 3 hardly affects the calculated PMF profiles (Supplementary Figure S2D), which indicates that the two parallel 1 μs simulations are adequate for the CCK2R-G5 systems. For the CCK2R-G5NH2 system, the PMF profiles and the time evolutions of two reaction coordinates calculated from each trajectory are much similar (Supplementary Figures S2, S3). Both reaction coordinates for two trajectories have stable fluctuations in the whole 1 μs simulation time, which indicates the initial conformational state was very steady. However, to further verify the minimal PMF region for the CCK2R-G5NH2 system, we performed another two parallel 700 ns simulations (termed trajectory 3 and 4) from the same initial conformational state with different randomized initial atomic velocities. We can find that the minimal PMF region and the time evolutions of two reaction coordinates are much similar for four parallel trajectories (Supplementary Figure S2E–H). Moreover, adding trajectories 3 and 4 also hardly affects the calculated PMF profiles (Supplementary Figure S2I), which indicates that the two parallel 1 μs simulations are also adequate for the CCK2R-G5NH2 systems.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1054575/fphar-13-1054575-HTML/image_m/fphar-13-1054575-g001.jpg)
FIGURE 1. 2D PMF profiles of CCK2R-G5 (A) and CCK2R-G5NH2 (B) along the reaction coordinates of distances between the N/C-terminus of G5/G5NH2 (Cα atoms of G13 and F17, respectively) and the geometrical center of TM domain of CCK2R; The initial site from the cryo-EM structure was marked as symbol X (coordinate: 22.1, 10.8). The PMF profile of each system was calculated by reweighting the combined two 1 μs GaMD simulations.
To see how the complex structures changed along with the simulation time, we depicted the time evolutions of RMSD values of simulation structures of CCK2R (based on TM backbone atoms) and G5/G5NH2 (based on backbone atoms) from their initial structures (aligned with TM backbone atoms, Supplementary Figure S4). We can find that structures of CCK2R and G5NH2 were relatively stable in the whole 1 μs simulations in each trajectory, while structures of G5 endured larger changes, especially for those in CCK2R-G5 trajectory 1. Those results also indicated that the initial conformational state for CCK2R-G5NH2 was in an energy-stable state but for CCK2R-G5 was in an energy-unstable state.
To make a detailed comparison between the binding modes and the underlying mechanisms of CCK2R in complexes with G5 and G5NH2, we extracted about 5,000 frames with PMF values equal to 0 kcal/mol evenly from the two parallel 1 μs simulations for the following analysis.
3.2 Conformational dynamics of CCK2R in complex with G5/G5NH2
PCA was performed (based on the Cα atoms of CCK2R) to investigate whether the primary conformational features of CCK2R in complexes with G5 and G5NH2 were the same. To make the dynamics of CCK2R conformations can be compared within a common subspace, the CCK2R conformations extracted from the CCK2R-G5 and CCK2R-G5NH2 complexes were combined and aligned with the cryo-EM structure based on the Cα atoms of TM domain first. The projections of CCK2R conformations onto the first two PCs were presented in Figure 2A. We can see that an obvious difference occurred at PC1 (the most important principal component of conformational dynamics), with negative and positive values for CCK2R in complexes with G5NH2 and G5, respectively. Different PC1 values indicate that the primary conformational features of CCK2R binding with G5 and G5NH2 are significantly different. However, CCK2R in cryo-EM structure (binding with gastrin-17) and CCK2R-G5NH2 simulation structures were located in the same region (Figure 2A), indicating that CCK2R binding with gastrin-17 and G5NH2 has similar primary conformational features. To clearly illustrate the primary structural differences between CCK2R binding with G5 and G5NH2, structural motions of CCK2R along with the PC1 were depicted in Supplementary Figure S5A. Form Supplementary Figure S5A, we can find that relative to binding with G5NH2, the upward parts of TM1, 2 and 5 of CCK2R binding with G5 moved closer to the center of the TM bundle. In contrast, other TMs moved away from the center of the TM bundle. Conformations with median PC1 and PC2 values were selected as the representative conformations and were depicted in Figure 2B. We can see that G5NH2 stayed in the TM domain of CCK2R while G5 moved out to the extracellular domain. The departure of G5 from the TM domain may be why the upward parts of TM1 and TM2 moved closer to the center of the TM bundle. Besides, conformations of the CCK2R-G5NH2 complex obtained in our GaMD simulations with minimal PMF values were much similar to the corresponding parts in the cryo-EM structure (Supplementary Figure S5B). The mean RMSDs of CCK2R (calculated on the TM backbone atoms) and G5NH2 (calculated on backbone atoms) from the cryo-EM structure (aligned with the TM backbone atoms) were 2.37 ± 0.39 Å and 1.70 ± 0.35 Å, respectively. At the same time, the minimal RMSD values were 1.51 Å and 1.25 for CCK2R and G5NH2, respectively.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1054575/fphar-13-1054575-HTML/image_m/fphar-13-1054575-g002.jpg)
FIGURE 2. (A) Projections of CCK2R conformations in the complexes of CCK2R-G5 (cyan dots) and CCK2R-G5NH2 (orange triangles) with PMF equal to 0 kcal/mol onto the common subspace defined by the first two PCs. The representative conformations in B are marked as symbol X, and the cryo-EM structure is marked as symbol Y. (B) Comparison of representative conformations of CCK2R-G5 (cyan) and CCK2R-G5NH2 (orange) complexes with PMF equal to 0 kcal/mol. CCK2R is shown as cartoons. G5 and G5NH2 are shown as cartoons, and their C-termini are shown as sticks. The directions of TM movement of CCK2R along with the PC1 from lowest to highest are indicated by arrows. (C) Linear mutual information (LMI) correlations between Cα atom pairs within CCK2R-G5 (upper triangle) and CCK2R-G5NH2 (lower triangle) complexes. Residues from 55 to 405 belong to CCK2R, and the rest 5 (at the end of the coordinate axis) belong to G5/G5NH2. (D) Differences in correlation coefficients between CCK2R-G5NH2 and CCK2R-G5. The correlation coefficient is color-coded (see color bars).
Correlated motions of residue pairs in the complexes were calculated by the Linear Mutual Information (LMI) method and depicted in Figure 2C. For clarity, the differences in correlation coefficients between CCK2R-G5NH2 and CCK2R-G5 were depicted in Figure 2D. LMI has the merit of no unwanted dependency on the relative orientation of the fluctuations which the Pearson coefficient suffers from (Grant et al., 2021). From Figure 2C, D, we can find that residue pairs have stronger correlated motions in the CCK2R-G5NH2 complex than in the CCK2R-G5 complex, which indicates that structural fluctuations induced by G5NH2 binding are much easier to transmit to the intracellular parts than by G5 binding.
3.3 Binding free energies of G5/G5NH2 to CCK2R
The total binding free energies and the contributions of each energy component of G5/G5NH2 to CCK2R were calculated by the MM-PBSA method (Miller et al., 2012) and listed in Table 1. By investigating the contributions of each energy component listed in Table 1, we found that the hydrophobic interactions (comprised of ΔGgas-vdw and ΔGsol-nonpolar) played dominant roles in the total binding free energies (ΔGbinding) both for G5 and G5NH2. The dominant role of the hydrophobic interactions in the total binding free energy is also found in other protein-protein complexes formation (Simões et al., 2017; Macalino et al., 2018). To further investigate why G5 moved away from the preferred binding sites for G5NH2, we performed MM-PBSA analysis on the CCK2R-G5 complexes with similar binding modes as CCK2R-G5NH2 complexes. The corresponding results are listed in Supplymentary Table S2. We can find that for the similar binding modes, the hydrophobic interaction for binding with G5 is 12.23 kcal/mol larger than binding with G5NH2. This difference was much more significant than the difference in electrostatic interactions (−3.70 kcal/mol, comprised of ΔGgas-ele and ΔGsol-polar), which indicates that the TM bundle is unfavorable for ligands binding with many net negative charges.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1054575/fphar-13-1054575-HTML/image_m/fphar-13-1054575-t001.jpg)
TABLE 1. MM-PBSA derived binding free energies (in unit of kcal/mol) of G5 and G5NH2 to CCK2R. MM-PBSA was calculated on the conformations with minimal PMF value.
3.4 Mechanisms of different binding modes of CCK2R with G5 and G5NH2
To shed light on the dominant interactions responsible for the different binding modes of CCK2R with G5 and G5NH2, we decomposed the total binding free energy into the per-residue level (Figure 3). Due to the negatively charged pentapeptide of G5 and G5NH2, the absolute contributions of all charged residues in CCK2R are larger than 1 kcal/mol and, positively charged residues made favorable contributions, negatively charged residues made unfavorable contributions, respectively. Besides, due to the one more negative charge in G5, charged residues in CCK2R had apparent larger contributions (favorable or unfavorable) for binding with G5 than G5NH2 (Figures 3A–C). Another difference is the largest favorable contributions made by charged residues in CCK2R are Arg208ECL2 and Arg3566.58 for binding with G5 and G5NH2, respectively, which also indicated different binding sites for G5 and G5NH2. Besides, there are three uncharged resides (Pro197ECL2, Trp209ECL2, Asn3536.55) in CCK2R that made unfavorable contributions larger than 1 kcal/mol for binding with G5 but none for binding with G5NH2 (Figures 3B–D), which indicated there existed large unfavorable interactions between CCK2R and G5. By inspecting the per-residue contribution in G5 and G5NH2 (Supplymentary Table S3), we found that large differences occurred at the C-terminal two residues of Asp and Phe. The contribution of Asp in G5NH2 is more than 20 kcal/mol larger than in G5, while the contribution of the non-amidated Phe in G5 is more than 150 kcal/mol larger than the amidated Phe in G5NH2.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1054575/fphar-13-1054575-HTML/image_m/fphar-13-1054575-g003.jpg)
FIGURE 3. (A–D) Free energy (in unit of kcal/mol) contributions of residues in CCK2R for binding with G5 (A,B) and G5NH2 (C,D). Charged residues (Asp, Glu, Arg and Lys) in CCK2R (A,C) located in the binding sites were labeled. Uncharged residues in CCK2R (B,D) with absolute values of binding free energy contributions larger than 1 kcal/mol were labeled, in which residues located in and outside binding sites were labeled in black and red, respectively. Dashed lines in A and C, B, and D are the reference lines with values of ±5 and ±1, respectively. Detailed binding poses of CCK2R-G5 (E) and CCK2R-G5NH2 (F). Residues of G5 and G5NH2 are shown in stick representation. Residues of CCK2R are shown in line representation, Salt bridges (R208ECL2-F17 in E, R3566.58-D16 and R2155.36-D16 in F) and hydrogen bonds (N1152.65-M15, and Q204ECL2-G13 in E, Y1894.60-D16 and H207ECL2-D16 in F) between CCK2R and G5/G5NH2 are depicted by yellow and red dotted lines, respectively.
By observing the distributions of charged residues in CCK2R and G5/G5NH2 complexes, we found that relative to the single charged residue of Asp in G5NH2, charged residues of Asp and Phe in G5 were much closer to the positively charged Arg cluster (208ECL2, 2135.33, 2155.35, 3566.58, 365ECL3) in CCK2R (Supplementary Figure S6). In particular, Arg208ECL2 changed to the flat state pointing to G5 from the erect state pointing to solvent. By further investigating the atomic interactions between CCK2R and G5/G5NH2 (Figures 3E, F), we found that Arg208ECL2 formed a salt bridge with the C-terminal carboxyl of Phe in G5 but no residues formed salt bridges with the sidechain carboxyl of Asp, while Arg3566.58 and Arg2155.35 formed salt bridges with the sidechain carboxyl of Asp in G5NH2. However, Arg3566.58 and 2155.35 were also located in the binding sites near Phe in G5 (Figure 3E), and the total contributions made by them two were larger for binding with G5 (−81.65 kcal/mol) than binding with G5NH2 (−73.35 kcal/mol). The strong electrostatic attractions between the Arg cluster in CCK2R and Phe and Asp in G5 made the C-terminus of G5 move out from the initial binding sites, which were favorable for G5NH2. Besides, Gln204ECL2 and Asn1152.65 formed hydrogen bonds with Gly and Met of G5, respectively, both of which made the top three largest favorable contributions among uncharged residues for CCK2R binding with G5. The time evolutions of distances between atomic pairs involved in the three key interactions between CCK2R and G5 also showed that those three key interactions were very stable in the last 600 ns simulation time in trajectory 1 (Supplementary Figure S7E), which indicated that they played key roles in the stabilize the conformation of G5. However, these binding modes counted against forming salt bridges or hydrogen bonds between Arg208ECL2, 2155.35, 3566.58, Tyr1894.60, and His207ECL2 (forming hydrogen bonds with Asp in CCK2R-G5NH2, see below) and Asp in G5 (Supplementary Figure S7).
G5NH2 includes the essential bioactive unit of gastrin activating CCK2R. Thus, the interaction mechanism between G5NH2 and CCK2R can help us understand the activation mechanism of CCK2R. Binding conformations with the lowest PMF values in GaMD simulations are much similar to the counterparts in the static cryo-EM conformation (Figure 1B and Supplementary Figure S5B). Arg3566.58 and Arg2155.35 formed salt bridges with Asp, making the largest two favorable contributions. 12 uncharged residues made favorable contributions with absolute values larger than 1 kcal/mol, 10 of them were located in the binding sites (Figures 3D,F) and another 2 (Val1383.36 and Leu2225.42) were in the vicinity of binding sites with occurrence rates of 0.7 and 0.2, respectively. Among the uncharged residues, Tyr1894.60 and His207ECL2 made the first two largest favorable contributions for binding and both of them formed hydrogen bonds with Asp in G5NH2. For 14 residues located in the binding sites (Figure 3F), only two made favorable contributions with absolute values less than 1 kcal/mol (-0.94 and -0.74 for Ser1313.29 and Leu367ECL3, respectively). Thus, 12 key residues (Met1343.32, Gly1353.33, Tyr1894.60, Gln204ECL2, Cys205ECL2, Val206ECL2, His207ECL2, Arg2155.35, Asn3536.55, Arg3566.58, His364ECL3 and Ile3727.35, see methods) for binding with G5NH2 are distinguished. Among them, Arg3566.58, Arg2155.35, Tyr1894.60, and His207ECL2 played dominant roles in constraining the position of Asp in G5NH2. Furthermore, the importance of Arg3566.58, Tyr1894.60, and His207ECL2 for CCK2R binding with gastrin was verified by the alanine scanning mutagenesis, which would completely abolish the binding of gastrin-17 (Zhang et al., 2021).
4 Conclusion
In this research, we performed extensive GaMD simulations to investigate the molecular determinants in the process of gastrin C-terminal pentapeptide amide end activating CCK2R. The C-terminal indispensable amide cap could significantly influence the binding modes of pentagastrin to CCK2R. Amidated pentagastrin has much similar binding modes as the counterpart in gastrin-17. Non-amidated pentagastrin moved its negatively charged C-terminus to interact with the positively charged Arg cluster (208ECL2, 2135.33, 2155.35, 3566.58, 365ECL3) distributed near the extracellular domain of CCK2R. Different binding modes induced different conformational dynamics of CCK2R. Residue pairs in CCK2R binding with amidated pentagastrin had stronger correlated motions than binding with non-amidated pentagastrin. 12 Key residues important for CCK2R-pentagastrin binding were identified. Among them, Arg3566.58, Arg2155.35, Tyr1894.60, and His207ECL2 formed strong salt bridges or hydrogen bonds with Asp in the amidated pentagastrin. Thus, they played dominant roles in constraining the position of Asp. In summary, our results explained the indispensable role of the C-terminal amide cap for gastrin’s bioactivity on a molecular basis. And the identified key interactions between the essential bioactive unit of gastrin and CCK2R would provide significant insights for developing CCK2R antagonists.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author contributions
KY, GW, and GZ contributed to the conception and design of the study. KY and XG performed the experiments. KY, HJ, GW, and GZ wrote the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This research was funded by the first batch of science and technology projects for the construction of innovative ecosystems by the National Supercomputing Center in Zhengzhou (grant No. 201400210500) and GHfund (grant No. 202202022169).
Acknowledgments
We are grateful to National Supercomputing Center in Zhengzhou for essential support.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.1054575/full#supplementary-material
References
Anders, J., Blüggel, M., Meyer, H. E., Kühne, R., ter Laak, A. M., Kojro, E., et al. (1999). Direct identification of the agonist binding site in the human brain CholecystokininB receptor. Biochemistry 38 (19), 6043–6055. doi:10.1021/bi990269z
Bhattarai, A., Wang, J., and Miao, Y. (2020). G-Protein-Coupled receptor–membrane interactions depend on the receptor activation state. J. Comput. Chem. 41 (5), 460–471. doi:10.1002/jcc.26082
Bläker, M., Ren, Y., Gordon, M. C., Hsu, J. E., Beinborn, M., and Kopin, A. S. (1998). Mutations within the cholecystokinin-B/gastrin receptor ligand ‘pocket’ interconvert the functions of nonpeptide agonists and antagonists. Mol. Pharmacol. 54 (5), 857–863. doi:10.1124/mol.54.5.857
Bläker, M., Ren, Y., Seshadri, L., McBride, E. W., Beinborn, M., and Kopin, A. S. (2000). CCK-B/Gastrin receptor transmembrane domain mutations selectively alter synthetic agonist efficacy without affecting the activity of endogenous peptides. Mol. Pharmacol. 58 (2), 399–406. doi:10.1124/mol.58.2.399
Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Ca. Cancer J. Clin. 68(6), 394–424. doi:10.3322/caac.21492
Calatayud, S., Alvarez, A., and Victor, V. M. (2010). Gastrin: An acid-releasing, proliferative and immunomodulatory peptide? Mini Rev. Med. Chem. 10(1), 8–19. doi:10.2174/138955710791112532
Cui, F. C., Yang, K. C., and Li, Y. Q. (2015). Investigate the binding of catechins to trypsin using docking and molecular dynamics simulation. Plos One 10 (5), e0125848. doi:10.1371/journal.pone.0125848
Dimaline, R., and Varro, A. (2014). Novel roles of gastrin. J. Physiol. 592 (14), 2951–2958. doi:10.1113/jphysiol.2014.272435
Ding, Y., Zhang, H., Liao, Y.-Y., Chen, L.-N., Ji, S.-Y., Qin, J., et al. (2022). Structural insights into human brain–gut peptide cholecystokinin receptors. Cell Discov. 8 (1), 55. doi:10.1038/s41421-022-00420-3
Dockray, G. J. (2004). Clinical endocrinology and metabolism Gastrin. Best. Pract. Res. Clin. Endocrinol. Metab. 18 (4), 555–568. doi:10.1016/j.beem.2004.07.003
Dolcetti, R., Re, V. D., and Canzonieri, V. (2018). Immunotherapy for gastric cancer: Time for a personalized approach? Int. J. Mol. Sci. 19 (6), 1602. doi:10.3390/ijms19061602
Eastman, P., Swails, J., Chodera, J. D., McGibbon, R. T., Zhao, Y., Beauchamp, K. A., et al. (2017). OpenMM 7: Rapid development of high performance algorithms for molecular dynamics. PLoS Comput. Biol. 13 (7), e1005659. doi:10.1371/journal.pcbi.1005659
Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103 (19), 8577–8593. doi:10.1063/1.470117
Feller, S. E., Zhang, Y., Pastor, R. W., and Brooks, B. R. (1995). Constant pressure molecular dynamics simulation: The Langevin piston method. J. Chem. Phys. 103(11), 4613–4621. doi:10.1063/1.470648
Foucaud, M., Archer-Lahlou, E., Marco, E., Tikhonova, I. G., Maigret, B., Escrieut, C., et al. (2008). Insights into the binding and activation sites of the receptors for cholecystokinin and gastrin. Regul. Pept. 145 (1-3), 17–23. doi:10.1016/j.regpep.2007.09.024
Fourmy, D., Gigoux, V., and Reubi, J. C. (2011). Gastrin in gastrointestinal diseases. Gastroenterology 141 (3), 814–818. e813. doi:10.1053/j.gastro.2011.07.006
Galés, C., Poirot, M., Taillefer, J., Maigret, B., Martinez, J., Moroder, L., et al. (2003). Identification of tyrosine 189 and asparagine 358 of the cholecystokinin 2 receptor in direct interaction with the crucial C-terminal amide of cholecystokinin by molecular modeling, site-directed mutagenesis, and structure/affinity studies. Mol. Pharmacol. 63 (5), 973–982. doi:10.1124/mol.63.5.973
Gilliam, A. D., and Watson, S. A. (2007). G17DT: An antigastrin immunogen for the treatment of gastrointestinal malignancy. Expert Opin. Biol. Ther. 7 (3), 397–404. doi:10.1517/14712598.7.3.397
Grant, B. J., Skjærven, L., and Yao, X.-Q. (2021). The Bio3D packages for structural bioinformatics. Protein Sci. 30 (1), 20–30. doi:10.1002/pro.3923
Hou, T., Wang, J., Li, Y., and Wang, W. (2011). Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations. J. Chem. Inf. Model. 51 (1), 69–82. doi:10.1021/ci100275a
Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., de Groot, B. L., et al. (2017). CHARMM36m: An improved force field for folded and intrinsically disordered proteins. Nat. Methods 14 (1), 71–73. doi:10.1038/nmeth.4067
Hussain, S. P. (2016). Pancreatic cancer: Current progress and future challenges. Int. J. Biol. Sci. 12 (3), 270–272. doi:10.7150/ijbs.14950
Irina, G. T., Esther, M., Elodie, L.-A., Ingrid, L., Magali, F., Bernard, M., et al. (2007). Validated ligand binding sites in CCK receptors. Next step: Computer-aided design of novel CCK ligands. Curr. Top. Med. Chem. 7 (12), 1243–1247. doi:10.2174/156802607780960528
Isberg, V., de Graaf, C., Bortolato, A., Cherezov, V., Katritch, V., Marshall, F. H., et al. (2015). Generic GPCR residue numbers - aligning topology maps while minding the gaps. Trends Pharmacol. Sci. 36 (1), 22–31. doi:10.1016/j.tips.2014.11.001
Jo, S., Kim, T., Iyer, V. G., and Im, W. (2008). CHARMM-GUI: A web-based graphical user interface for CHARMM. J. Comput. Chem. 29 (11), 1859–1865. doi:10.1002/jcc.20945
Lange, O. F., and Grubmüller, H. (2006). Generalized correlation for biomolecular dynamics. Proteins 62 (4), 1053–1061. doi:10.1002/prot.20784
Langer, I., Tikhonova, I. G., Travers, M.-A., Archer-Lahlou, E., Escrieut, C., Maigret, B., et al. (2005). Evidence that interspecies polymorphism in the human and rat cholecystokinin receptor-2 affects structure of the binding site for the endogenous agonist cholecystokinin. J. Biol. Chem. 280 (23), 22198–22204. doi:10.1074/jbc.M501786200
Lindorff-Larsen, K., Piana, S., Dror, R. O., and Shaw, D. E. (2011). How fast-folding proteins fold. Science 334 (6055), 517–520. doi:10.1126/science.1208351
Lipiński, P. F. J., Garnuszek, P., Maurin, M., Stoll, R., Metzler-Nolte, N., Wodyński, A., et al. (2018). Structural studies on radiopharmaceutical DOTA-minigastrin analogue (CP04) complexes and their interaction with CCK2 receptor. EJNMMI Res. 8 (1), 33. doi:10.1186/s13550-018-0387-3
Macalino, S. J. Y., Basith, S., Clavio, N. A. B., Chang, H., Kang, S., and Choi, S. (2018). Evolution of in silico strategies for protein-protein interaction drug discovery. Molecules 23 (8), 1963. doi:10.3390/molecules23081963
Maddalo, G., Spolverato, Y., Rugge, M., and Farinati, F. (2014). Gastrin: From pathophysiology to cancer prevention and treatment. Eur. J. Cancer Prev. 23 (4), 258–263. doi:10.1097/cej.0000000000000008
Magnan, R., Escrieut, C., Gigoux, V., De, K., Clerc, P., Niu, F., et al. (2013). Distinct CCK-2 receptor conformations associated with beta-arrestin-2 recruitment or phospholipase-C activation revealed by a biased antagonist. J. Am. Chem. Soc. 135 (7), 2560–2573. doi:10.1021/ja308784w
Marco, E., Foucaud, M., Langer, I., Escrieut, C., Tikhonova, I. G., and Fourmy, D. (2007). Mechanism of activation of a G protein-coupled receptor, the human cholecystokinin-2 receptor. J. Biol. Chem. 282 (39), 28779–28790. doi:10.1074/jbc.M700349200
Miao, Y., and McCammon, J. A. (2016). G-Protein coupled receptors: Advances in simulation and drug discovery. Curr. Opin. Struct. Biol. 41, 83–89. doi:10.1016/j.sbi.2016.06.008
Miao, Y., and McCammon, J. A. (2017). “Gaussian accelerated molecular dynamics: Theory, implementation, and applications,”. Editor D. A. Dixon (Elsevier) 13, 231–278.Annu. Rep. Comput. Chem. doi:10.1016/bs.arcc.2017.06.005
Miao, Y., Sinko, W., Pierce, L., Bucher, D., Walker, R. C., and McCammon, J. A. (2014). Improved reweighting of accelerated molecular dynamics simulations for free energy calculation. J. Chem. Theory Comput. 10 (7), 2677–2689. doi:10.1021/ct500090q
Miller, B. R., McGee, T. D., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012). MMPBSA.py: An efficient program for end-state free energy calculations. J. Chem. Theory Comput. 8 (9), 3314–3321. doi:10.1021/ct300418h
Mori, T., Miyashita, N., Im, W., Feig, M., and Sugita, Y. (2016). Molecular dynamics simulations of biological membranes and membrane proteins using enhanced conformational sampling algorithms. Biochim. Biophys. Acta 1858 (7), 1635–1651. Part B. doi:10.1016/j.bbamem.2015.12.032
Palermo, G., Miao, Y., Walker, R. C., Jinek, M., and McCammon, J. A. (2017). CRISPR-Cas9 conformational activation as elucidated from enhanced molecular simulations. Proc. Natl. Acad. Sci. U. S. A. 114 (28), 7260–7265. doi:10.1073/pnas.1707645114
Pang, Y. T., Miao, Y., Wang, Y., and McCammon, J. A. (2017). Gaussian accelerated molecular dynamics in NAMD. J. Chem. Theory Comput. 13 (1), 9–19. doi:10.1021/acs.jctc.6b00931
Pannequin, J., Barnham, K. J., Hollande, F., Shulkes, A., Norton, R. S., and Baldwin, G. S. (2002). Ferric ions are essential for the biological activity of the hormone glycine-extended gastrin. J. Biol. Chem. 277 (50), 48602–48609. doi:10.1074/jbc.M208440200
Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., et al. (2005). Scalable molecular dynamics with NAMD. J. Comput. Chem. 26 (16), 1781–1802. doi:10.1002/jcc.20289
Rehfeld, J. F. (2017). Cholecystokinin-from local gut hormone to ubiquitous messenger. Front. Endocrinol. 8, 47. doi:10.3389/fendo.2017.00047
Ritler, A., Shoshan, M., Deupi, X., Wilhelm, P., Schibli, R., Wennemers, H., et al. (2019). Elucidating the structure-activity relationship of the pentaglutamic acid sequence of minigastrin with cholecystokinin receptor subtype 2. Bioconjug. Chem. 30, 657–666. doi:10.1021/acs.bioconjchem.8b00849
Seeber, M., Felline, A., Raimondi, F., Muff, S., Friedman, R., Rao, F., et al. (2011). Wordom: A user-friendly program for the analysis of molecular structures, trajectories, and free energy surfaces. J. Comput. Chem. 32 (6), 1183–1194. doi:10.1002/jcc.21688
Silvente-Poirot, S., Escrieut, C., Galès, C., Fehrentz, J.-A., Escherich, A., Wank, S. A., et al. (1999). Evidence for a direct interaction between the penultimate aspartic acid of cholecystokinin and histidine 207, located in the second extracellular loop of the cholecystokinin B receptor. J. Biol. Chem. 274 (33), 23191–23197. doi:10.1074/jbc.274.33.23191
Silvente-Poirot, S., Escrieut, C., and Wank, S. A. (1998). Role of the extracellular domains of the cholecystokinin receptor in agonist binding. Mol. Pharmacol. 54 (2), 364–371. doi:10.1124/mol.54.2.364
Silvente-Poirot, S., and Wank, S. A. (1996). A segment of five amino acids in the second extracellular loop of the cholecystokinin-B receptor is essential for selectivity of the peptide agonist gastrin. J. Biol. Chem. 271 (25), 14698–14706. doi:10.1074/jbc.271.25.14698
Simões, I. C. M., Costa, I. P. D., Coimbra, J. T. S., Ramos, M. J., and Fernandes, P. A. (2017). New parameters for higher accuracy in the computation of binding free energy differences upon alanine scanning mutagenesis on protein–protein interfaces. J. Chem. Inf. Model. 57 (1), 60–72. doi:10.1021/acs.jcim.6b00378
Sitarz, R., Skierucha, M., Mielko, J., Offerhaus, G. J. A., Maciejewski, R., and Polkowski, W. P. (2018). Gastric cancer: Epidemiology, prevention, classification, and treatment. Cancer Manag. Res. 10, 239–248. doi:10.2147/cmar.s149619
Smith, J. P., Fonkoua, L. K., and Moody, T. W. (2016). The role of gastrin and CCK receptors in pancreatic cancer and other malignancies. Int. J. Biol. Sci. 12 (3), 283–291. doi:10.7150/ijbs.14952
Song, J., Ren, H., Li, Y., Xu, J., Kong, H., Tong, W., et al. (2013). rG17PE38, a novel immunotoxin target to gastric cancer with overexpressed CCK-2R. J. Drug Target. 21 (4), 375–382. doi:10.3109/1061186x.2012.757770
Stone, S. R., Mierke, D. F., and Jackson, G. E. (2007). Evidence for a C-terminal structural motif in gastrin and its bioactive fragments in membrane mimetic media. Peptides 28 (8), 1561–1571. doi:10.1016/j.peptides.2007.07.009
Wang, C., Greene, D. A., Xiao, L., Qi, R., and Luo, R. (2018). Recent developments and applications of the MMPBSA method. Front. Mol. Biosci. 4 (87), 87. doi:10.3389/fmolb.2017.00087
Weinbach, Y., and Elber, R. (2005). Revisiting and parallelizing SHAKE. J. Comput. Phys. 209 (1), 193–206. doi:10.1016/j.jcp.2005.03.015
Willard, M. D., Lajiness, M. E., Wulur, I. H., Feng, B., Swearingen, M. L., Uhlik, M. T., et al. (2012). Somatic mutations in CCK2R alter receptor activity that promote oncogenic phenotypes. Mol. Cancer Res. 10 (6), 739–749. doi:10.1158/1541-7786.mcr-11-0483
Wu, E. L., Cheng, X., Jo, S., Rui, H., Song, K. C., Dávila-Contreras, E. M., et al. (2014). CHARMM-GUI Membrane Builder toward realistic biological membrane simulations. J. Comput. Chem. 35 (27), 1997–2004. doi:10.1002/jcc.23702
Yang, K.-C., Cui, F.-C., Shi, C., Chen, W.-D., and Li, Y.-Q. (2019). Evolution of conformation and dynamics of solvents in hydration shell along the urea-induced unfolding of ubiquitin. Chin. J. Polym. Sci. 37 (7), 708–718. doi:10.1007/s10118-019-2238-8
Yang, K., Cui, F., and Li, Y. (2018). Distribution and dynamics of water and urea in hydration shell of ribonuclease sa:A molecular dynamics simulation study. Chin. J. Appl. Chem. 35 (10), 1243–1248. doi:10.11944/j.issn.1000-0518.2018.10.170342
Yang, K., Różycki, B., Cui, F., Shi, C., Chen, W., and Li, Y. (2016). Sampling enrichment toward target structures using hybrid molecular dynamics-Monte Carlo simulations. PLoS ONE 11 (5), e0156043. doi:10.1371/journal.pone.0156043
Zhang, X., Cui, F., Chen, H., Zhang, T., Yang, K., Wang, Y., et al. (2018). Dissecting the innate immune recognition of opioid inactive isomer (+)-Naltrexone derived toll-like receptor 4 (TLR4) antagonists. J. Chem. Inf. Model. 58 (4), 816–825. doi:10.1021/acs.jcim.7b00717
Keywords: gastrin, cholecystokinin 2 receptor, interaction mechanism, GAMD, PMF, MM-PBSA analysis
Citation: Yang K, Jin H, Gao X, Wang G-C and Zhang G-Q (2023) Elucidating the molecular determinants in the process of gastrin C-terminal pentapeptide amide end activating cholecystokinin 2 receptor by Gaussian accelerated molecular dynamics simulations. Front. Pharmacol. 13:1054575. doi: 10.3389/fphar.2022.1054575
Received: 27 September 2022; Accepted: 02 December 2022;
Published: 23 January 2023.
Edited by:
Jinan Wang, University of Kansas, United StatesReviewed by:
David Bickel, Vrije University Brussels, BelgiumMarcel Bermudez, University of Münster, Germany
Copyright © 2023 Yang, Jin, Gao, Wang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kecheng Yang, yangkch@zzu.edu.cn