Molecular Dynamics and Bioactivity of a Novel Mutated Human Parathyroid Hormone

Purpose : To design and evaluate a novel human parathyroid hormone (hPTH) analog. Methods : Mutation stability prediction was processed on hPTH, docked the mutant hPTH with its receptor, and then proceeded with molecular dynamics using Discovery Studio 3.5 software package for the complex. The  bioactivity of the hPTH on the expression levels of the Rps27, RANKL and OPG genes were assessed in UAMS-32P cells. Results : A three-site mutant, hPTH (R25Q:K26E:K27L), was obtained and MD trajectory analysis showed a decrease in the root mean square deviation by 51.95%, in the radius of gyration by 3.57%, and in the interaction energy by 10%, compared with the wild-type hPTH. Bioactivity assessment demonstrated that this mutant stimulated the ratio of RANKL/OPG 30-fold higher than the wild-type. Conclusion : We successfully designed a new hPTH mutant with a stable conformation and high bioactivity, and this may be useful for elucidating ligand-receptor recognition mechanism and discovering novel hPTH analogs. Keywords : Parathyroid hormone, Mutation prediction, Molecular dynamics, RANKL/OPG, UAMS-32P cell


INTRODUCTION
Parathyroid hormone (PTH) is an 84-amino-acid polypeptide that regulates extracellular calcium homeostasis and bone remodeling by direct action on the kidneys and bones [1]. PTH exerts its biological effects by binding to and activating the PTH receptor (PTH1R), which is highly expressed in PTH target tissues [2]. As a member of family B G protein-coupled receptor (GPCR), the PTH1R contains an N-terminal extracellular domain (ECD) with three conserved disulfide bonds and a C-terminal domain with seven transmembrane helices [3]. The N-terminal 34-residue region of PTH (1-34) is capable of fully activating PTH1Rs to the same degree as PTH . Daily administration of either versions of PTH enhances bone formation, increases bone density, and decreases fracture risk in postmenopausal women [4,5].
As a therapeutic peptide, human PTH (1-34) (hPTH) has natural limitations, such as low oral bioavailability, short half-life (only 5 min in serum) [6], high conformational flexibility (NMR studies revealed that hPTH (1-34) has a U-shaped tertiary structure [7], whereas its X-ray crystal structure is a single continuous helix [8]), and high synthesis and production costs (Forteo, a recombinant hPTH of Eli Lilly, costs US $16.70 per day and US $6,100 for a year of optimal therapy), which is a great therapeutic challenge to health-care systems [9]. Production of hPTH with a more stable structure and a higher PTH1R binding affinity is one way to overcome these shortcomings.
Calculate Mutation Energy (Stability) protocol evaluates the effect of single-point mutations on protein stability. It performs amino-acid scanning mutagenesis on a set of selected amino-acid residues by mutating each of them to one or more specified amino-acid types. The energy effect of each mutation on the protein stability (mutation energy, ΔΔGmut) is calculated as the difference of the free energy of folding between the mutated structure and the wild type protein (Eq 1). ΔΔGmut = ΔΔGfold (mutant) -ΔΔGfold (wild type) ….. (1) where the folding free energy, ΔΔGfold, is defined as free energy difference between the folded and unfolded states of the protein. The unfolded state is modeled as a relaxed peptide in extended conformation with the mutated residue in the center. This is used to account for the van der Waals and short-range electrostatic interactions between the mutated residue and the rest of the protein. In addition, a Gaussian random chain model is used to estimate the long range electrostatic interaction energy of the unfolded state. All energy terms are calculated by CHARMm and the electrostatic energy is calculated using a Generalized Born implicit solvent model.

Molecular dynamics (MD) simulations
Complex of the wild-type WT-hPTH (15-34) with PTH1R and the mutant MT-hPTH (15-34) with PTH1R were used as MD simulation starting points. A full-atom system of PTH1R/WT-hPTH (15-34) in a solution box was constructed using the Discovery Studio 3.5 solvation protocol with default parameters. All ionizable residues were assigned to their respective charge status corresponding to pH 7.5. The complex was then solvated in a water box by the Explicit Periodic Boundary solvation model with a minimum distance of 7.0 Å from the boundary. The cell shape was orthorhombic with 116.686×88.7512×65.3422 Å 3 dimensions. Similarly, we also constructed an all-atom PTH1R/MT-hPTH (15-34) solution box system with 116.837×88.7548×65.5112 Å 3 cell dimensions.
The system was subjected to the CHARMm force-field and relaxed by energy minimization (2000 steps of steepest descent and 2000 steps of conjugated gradient) with the harmonic restraint for the protein. The system was slowly driven from an initial temperature of 50 K to the target temperature of 300 K for 10 ps and equilibration simulations were run for 10 ps. The MD simulations (production) were performed for 10 ns with the NPT system at a constant temperature of 300 K and the results were saved at a frequency of 1000 steps. All the other parameters were set as defaults.

MD trajectory analysis
The MD trajectory was determined for structural properties, root mean-square deviation (RMSD), potential energy, and radius of gyration (Rg) using the Discovery Studio 3.5 analyze trajectory protocol. The interaction energy between hPTH and the PTH1R was then calculated using the calculate interaction energy protocol. All parameters were set at defaults.

Peptide synthesis
The WT-hPTH (1-34) was obtained from Sigma Corp. (Santa Clara, CA). The MT-hPTH (1-34) was prepared by the solid-phase procedure and purified by gel filtration and sequential HPLC by GLS Biochem Co, Ltd. (Shanghai, China).

Cell culture and bioactivity evaluation
Confluent UAMS-32P cells stably expressing PTH1R were used to evaluate the bioactivity of WT-hPTH (1-34) and synthesized MT-hPTH (1-34). The effects of the hPTHs on the expression levels of the Rps27, RANKL and OPG genes were assessed using RT-PCR.
The UAMS-32P cell line was maintained in αminimal essential medium (Invitrogen, Carlsbad, CA) supplemented with 10% fetal bovine serum (FBS, HyClone, Logan, UT) and 1% each of penicillin, streptomycin, and glutamine. The vehicle was set as the negative control, along with WT-hPTH (1-34) as positive control, and the synthesized MT-hPTH (1-34) as the sample of interest. The UAMS-3P cells were treated with each agent at a concentration of 1×10 -7 mol/L and total RNA was extracted from the cells after 24 h.

Statistical analysis
The software suite of life science molecular design solutions Discovery Studio (version 3.5, Accelrys Inc., San Diego, CA) was used for protein design and molecular dynamics simulation. The analysis of the MD data was made using Origin Pro (version 9.0). Mean ± standard error of mean (SEM) of the data were computed. The experimental data of bioactivity in Fig 3 were analyzed with GraphPad Prism (version 5.0). The level of significance was p < 0.01.

RESULTS hPTH mutation prediction
We carried out full mutation scanning on Leu15 to Phe34 and obtained 380 single mutants, 1075 double mutants, and 637 triple mutants, then sorted them by their mutation energy. The two mutants with the lowest mutation energy of each type are listed in Table 1. Mutation energy prediction showed that the triple mutation had a much lower energy than the other two types. We then selected 5 mutants with the lowest mutation energy from each mutation type for further study (data not shown). In this study, the hPTH with the triple mutation of Arg25Gln, Lys26Glu, Lys27Leu was chosen for analysis.

MD simulation trajectory
The structural changes and the binding between hPTH and PTH1R were investigated by viewing the system potential energy, Rg, RMSD, and interaction energy (Fig 1). Statistical analysis was carried out for the last 5 ns of MD simulations to correlate the MD data with peptide characteristics ( Table 2).
Hydrophobicity analysis showed that hydrophobic residues of MT-hPTH (15-34) lay along the amphipathic helix and faced a hydrophobic groove formed in the receptor interface. With the aid of hydrophobic force, Phe34 of the mutant was linked to Lys158. But Val31 and Phe34 of the wild-type did not establish contact with the receptor through hydrophobic interaction since they were in the turn or coil region of the peptide (Fig 2b & d).
To investigate the structural fluctuation, the RMSF of each residue of WT-hPTH (15-34) and MT-hPTH (15-34) was calculated (Fig 2e). The results showed that the RMSF values of the residues of the mutant were low and relatively constant, while those of the wild-type ranged from 0.2 to 1.1 nm. The flexible residues in the 29-34 region were located at the C-terminal of the wild-type. We demonstrated that the Arg25Gln, Lys26Glu, and Lys27Leu mutations strengthened the stability of the C-terminal region of hPTH.

MT-hPTH (1-34) bioactivity
PTH (1-34) stimulates the transcription of the RANKL gene and inhibits OPG gene transcription, increasing the RANKL/OPG ratio. The ratio of RANKL/OPG expression by osteoblasts is reported to be a key determinant of osteoclastogenic activity [14].

DISCUSSION
Leu24 and Leu28 of hPTH (1-34) are located at the center of the hydrophobic interface [8]. Sitedirected mutagenesis studies on the C-terminal region of hPTH (1-34) have suggested that Leu24 and Leu28 are intolerant to mutation. When Leu24 and Leu28 are substituted by Glu, the receptor binding affinity decreases by 4000and 1600-fold, respectively [15]. In our mutation prediction, the mutation energy of Leu24 and Leu28 was high and the mutation effect was unstable (data not shown), which is largely in agreement with the previous study. Arg25, Lys26, and Lys27 located between the intolerant residues Leu24 and Leu28 are mutation-tolerant [16].
In addition, we found that the Lys26-Glu mutation led to the formation of two new H-bond pairs. Lys27 is hydrophilic polar amino-acid, whereas its mutant residue Leu is hydrophobic and nonpolar, so this mutation increased the hydrophobic area between hPTH and its receptor (Fig 2b & d). It is likely that the increase of H-bond and hydrophobic force resulted in the enhancement of ligand-receptor interaction energy.
So far, mutagenesis studies of PTH have been carried out in PTH-specific binding sites and their interaction with its receptor [11,17]; however there are no reports on the prediction of hPTH multi-site mutations. We performed mutation prediction and obtained a triple hPTH mutant, hPTH (R25Q:K26E:K27L), which evidently increased the RANKL/OPG ratio by 30 times, but its mechanism of action remains unclear. These findings therefore present a new conceptual starting-point for unraveling the ligand-receptor recognition mechanism and consequently provide a guide to designing novel hPTH analogs.

CONCLUSION
PTH is a key regulator of calcium homeostasis and a major mediator of bone remodeling. Due to its natural limitations, we carried out mutation prediction and bioactivity evaluation to get a mutant with higher bioactivity and more structural stability.
The ligand-receptor recognition mechanism was also investigated. As a result, we successfully designed a new hPTH mutant with a stable conformation and high bioactivity, which may be useful for elucidating the ligandreceptor recognition mechanism and discovering novel hPTH analogs.