Original Research ARTICLE
Exploring the Interaction Mechanism Between Cyclopeptide DC3 and Androgen Receptor Using Molecular Dynamics Simulations and Free Energy Calculations
- 1School of Pharmacy, Lanzhou University, Lanzhou, China
- 2College of Chemistry and Chemical Engineering, Lanzhou University, Lanzhou, China
Androgen receptor (AR) is a key target in the discovery of anti-PCa (Prostate Cancer) drugs. Recently, a novel cyclopeptide Diffusa Cyclotide-3 (DC3), isolated from Hedyotisdiffusa, has been experimentally demonstrated to inhibit the survival and growth of LNCap cells, which typically express T877A-mutated AR, the most frequently detected point mutation of AR in castration-resistant prostate cancer (CRPC). But the interaction mechanism between DC3 and AR is not clear. Here in this study we aim to explore the possible binding mode of DC3 to T877A-mutated AR from molecular perspective. Firstly, homology modeling was employed to construct the three-dimensional structure of the cyclopeptide DC3 using 2kux.1.A as the template. Then molecular docking, molecular dynamics (MD) simulations, and molecular mechanics/generalized Born surface area (MM-GBSA) methods were performed to determine the bind site and explore the detailed interaction mechanism of DC3-AR complex. The obtained results suggested that the site formed by H11, loop888-893, and H12 (site 2) was the most possible position of DC3 binding to AR. Besides, hydrogen bonds, hydrophobic, and electrostatic interactions play dominant roles in the recognition and combination of DC3-AR complex. The essential residues dominant in each interaction were specifically revealed. This work facilitates our understanding of the interaction mechanism of DC3 binding to AR at the molecular level and contributes to the rational cyclopeptide drug design for prostate cancer.
Prostate cancer (PCa) has become the second frequently diagnosed cancer in men throughout the world (American Cancer Society, 2015). Prostate, lung and bronchus, colorectal cancers accounts for about 44% of all cancer cases in men, with PCa alone accounting for 1 in 5 new diagnoses (Siegel et al., 2016). PCa is especially common in economically developed countries and regions like Northern and Western Europe, Northern America, and Oceania (American Cancer Society, 2015). In America, prostate cancer is the most common cancer and was predicted as the leading cause of male cancer-related death over the next decade (Siegel et al., 2016). In those less developed countries, the incidence rate of prostate cancer is increasing with stable or increasing mortality trend in recent years (Center et al., 2012).
Androgen receptor (AR) (NR3C4, nuclear receptor subfamily 3, group C, gene 4), a member of steroid hormone group of nuclear receptor superfamily, plays an essential role in the development and proliferation of prostate cancer (Tsai and O'Malley, 1994; Mangelsdorf et al., 1995; Nuclear Receptors Nomenclature Committee, 1999). The survival and growth of PCa cells are dependent on the androgenic stimulation through AR. Firstly, 5α-dihydrotestosterone (DHT) binds to AR to promote the association of AR co-regulators. Then the activated AR migrates into nucleus and regulates the expression of target genes in prostate cells (Heinlein and Chang, 2004).
Clinically, PCa is commonly treated by AR pathway perturbation, such as androgen suppression via surgical or chemical castration [gonadotropin-releasing hormone (GnRH) analogs] means (Palmbos and Hussain, 2013). AR antagonist drugs, such as flutamide, nilutamide, bicalutamide, and enzalutamide, take effects by suppressing the action of androgens via competing for AR binding sites (Yamamoto et al., 2012). These androgen blockade therapies are initially effective, however, a considerable population of patients ultimately develop as castration-resistant prostate cancer (CRPC) after prolonged use of an AR antagonist (Schröder, 2008; Yamaoka et al., 2010). AR mutation is one of the leading causes of antiandrogens resistance (Tan et al., 2015). These mutated ARs bind to other steroid hormones and induce the activation of AR transcriptional activity in response to antiandrogens, which results in the PCa growth (Tan et al., 2015). In this case, it shows far-reaching significance to seek and explore novel anti-CRPC drugs targeting gene mutational AR.
DC3 (Diffusa Cyclotide-3) is a novel cyclopeptide isolated from the traditional Chinese Medicine (TCM) Hedyotisdiffusa (Hu et al., 2015), which has been widely used for the treatment of various cancers and tumors, including prostate cancer, in China with a long history (Lin et al., 2010, 2011; Liu et al., 2010; Lee et al., 2011). It has been experimentally detected that DC3 expresses potent cytotoxicity against LNCaP cells and inhibits the cell migration and invasion. Besides, it can significantly inhibit the development of tumor in weight and size in the mouse xenograft model. All these findings lead to the conclusion that DC3 has evident anti-PCa effects both in vitro and vivo (Hu et al., 2015). Moreover, in the DC3 sensitivity experiments on three types of human prostate cancer cells, androgen dependent LNCaP cell lines showed obvious higher sensitivity to DC3 comparing to androgen independent PC3 and DUl45 cell lines. Besides, LNCap cell lines typically express T877A-mutated AR, which is the most frequently detected point mutation in CRPC (Veldscholte et al., 1992; Zhou et al., 2010; Yamada et al., 2013). All these evidences suggest that DC3 is a potential candidate binding to T877A AR.
However, the interaction mechanisms between DC3 and T877A-mutated AR are not clear. Therefore, it will be constructive and profoundly significant to launch the mechanism-relevant research. Fortunately, the amino acid sequence of DC3 has been experimentally determined (Hu et al., 2015), which makes it possible to explore the interaction mechanism at the molecular level. Referencing to the published papers, the current research could be consist of three parts: (1) constructing the three-dimensional structure of the cyclopeptide (Jitonnom et al., 2012), in our case is DC3; (2) determining the binding site and binding pose of cyclopeptide to protein (Punkvang et al., 2015); (3) investigating the detailed interaction mechanism between cyclopeptide and protein (Liu et al., 2015; Hitzenberger et al., 2017).
In our study, first of all, homology modeling technology was conducted to construct the three-dimensional structure of cyclopeptide DC3 based on its amino acid sequence. Then molecular docking, all-atom molecular dynamics (MD) simulations and molecular mechanics/generalized Born surface area (MM/GBSA) methods and various MD trajectory analysis methods were combined to explore the most possible binding site of DC3 to AR, investigate the key residues dominant in the binding process, and elucidate the detailed interaction mechanism. The results are expected to reveal the interaction mechanism of DC3-AR complex, promote the development of DC3 and correlative cyclopeptide AR antagonist, which will contribute to the rational drug design for prostate cancer.
Homology Modeling of Cyclopeptide DC3
Homology modeling is a common technique to construct three-dimensional structure from amino acid sequence using homologous proteins with known structure as templates (Topham et al., 1990; Bordoli et al., 2009; Wang Z. et al., 2015). As amino acid sequence of DC3 was confirmed by Edman degradation and gene cloning (Hu et al., 2015), homology modeling was adopted here to build the 3D structure of DC3 using SWISS-MODEL (Arnold et al., 2006; Guex et al., 2009; Kiefer et al., 2009; Biasini et al., 2014; Bienert et al., 2016). Here, the SWISS-MODEL Template Library (SMTL) is searched both with BLAST and HHblits to identify templates and target-template alignments (Arnold et al., 2006). Then the template was selected based on various criteria such as sequence similarity, sequence identity, coverage, the global quality estimation score (GMQE) and so on.
Molecular Docking Analysis of DC3 to AR
Molecular docking (Benkert et al., 2011; Meng et al., 2011; Yuriev et al., 2015) was used to analyze the possible binding site and preferred orientation of DC3 into androgen receptor by simulating combining conformation and computing binding affinity. Here the crystal structure coordinates of the T877A-mutated AR LBD was obtained from the RCSB Protein Data Bank (http://www.rcsb.org/pdb; PDB ID: 4OHA). The missing loop regions were refined by Discovery Studio 2.5. (Accelrys Inc. CA, 2009). Molecular docking process was carried out by using ZDOCK module. The rigid-body protein–protein docking program ZDOCK uses the Fast Fourier Transform algorithm to enable an efficient global docking search on a 3D grid, and utilizes a combination of shape complementarity, electrostatics and statistical potential terms for scoring. Finally, two simple scoring functions–ZRank Score and ZDock Score, pose amount of each cluster, and the rationality of binding mode were taken into consideration to evaluate the docking results.
Molecular Dynamics Simulations
Molecular dynamics (MD) simulations were operated through Amber12 package (Case et al., 2012). All the simulations are under the circumstance of ff99SB force field (Hornak et al., 2006) and periodic boundary condition. Firstly, six chloride counterions were added to each system to maintain the electro-neutrality. Then all studied systems were, respectively, immersed into a cubic box of TIP3P (van der Spoel and van Maaren, 2006) water with edge of the box at least 10Å distant from the complex. Energy minimization was carried out in three stages with different harmonic restraint: all atoms constrained by 5.0 kcal·mol−1·Å−2, only receptor backbone atoms constrained by 3.0 kcal·mol−1·Å−2 and without any restraint. Each minimization was executed for 5,000 steps, in which the first 2,500 steps were calculated by the steepest descent method while the subsequent 2,500 steps were executed by conjugated gradient method. These systems were heated up to 310.0 K in the NVT ensemble for 100 ps with the receptor backbone atoms constrained by 5.0 kcal·mol−1·Å−2. And then, a total of 1.5 ns equilibration of each system was performed in NPT ensemble, where the former 800 ps were divided into four stages and the restraints applied to these stages were in a descending order (4.0, 3.0, 2.0, 1.0 kcal·mol−1·Å−2 respectively), the latter 700 ps were carried out without any restraint. Minimization and heat, as well as equilibration were executed in the Sander program. Finally, a 150 ns production of MD simulations of each system was performed in the PMEMD program at 310.0 K, 1 atm in the NPT ensemble without any restraint. The Langevin thermostat was used to control the temperature and the Berendsen barostat was used for constant pressure simulation. The time step was set as 2fs, and the coordinates of trajectories were recorded every 2 ps. During this simulation, the SHAKE algorithm (Ryckaert et al., 1977) was employed to constrain the bond lengths involving hydrogen, the Particle Mesh Ewald (PME) (Darden et al., 1993; Fischer et al., 2015) was adopted to calculate of electrostatic interaction with a 10Å non-bonded cutoff.
Free Energy Calculations
To investigate the interaction of DC3-AR systems from the energetic perspective, the binding free energy calculations based on the trajectories of MD simulations were performed by MM-GBSA method (Hou et al., 2011a,b; Xu et al., 2013; Sun et al., 2014a,b; Chen et al., 2016). The binding free energies ΔGbind was calculated as following equation:
where ΔH represents enthalpy contribution, which is composed of enthalpy changes in gas-phase (Egas) and solvent-phase (ΔGsol). –TΔS represents entropy contribution. Entropic calculation is time-consuming, and its value will fluctuate if a small quantity of snapshots were adopted (Hou et al., 2011a; Wang Q. et al., 2015). In this study it was omitted. Egas was considered as the sum of internal interaction (ΔEint) from bonds, angles, and torsions, van der Waals (ΔEvdw) and electrostatic energies (ΔEele) as follow:
ΔGsol can be decomposed into the polar and nonpolar contributions as follow:
Here, ΔGGB represents the polar solvation contribution, which is calculated by solving GB equation (Kollman et al., 2000; Onufriev et al., 2004). ΔGSA, estimated by the solvent accessible surface area, represents the nonpolar solvation contribution.
To further explore the detailed interaction information of DC3-AR complex, free energy decomposition was performed by using MM-GBSA method to identify the key residues responsible for binding energy. The contribution of each residue was calculated without considering the contribution of entropies. The contribution is defined as the sum of van der Waals contribution (ΔEvdW), electrostatic contribution (ΔEele), polar solvation contribution (ΔGGB), and nonpolar solvation contribution (ΔGSA):
Snapshots, used for both binding free energy and free energy decomposition calculations, were extracted from the last 50 ns of MD trajectories at intervals of 2 ps.
MD Trajectory Analysis
Hydrogen Bond Analysis
The number of formed hydrogen bonds vs. simulation time was calculated to detect the system stability during the process of simulation. Here, the hydrogen bond criteria was set as the distance of acceptor-donor <0.35 nm and the angle >120° (Fu et al., 2013). 0.35 nm is a common choice of hydrogen bond distance in literature (Liu et al., 2014; Wang Q. et al., 2015). The frames adopted for this calculation were extracted from the whole 150 ns MD trajectories at intervals of 2 ps. Besides, in order to determine exactly how hydrogen bonds play dominant roles in maintaining system stability in the last 50 ns MD simulations, the occupations of hydrogen bonds formed in this period were calculated as following equation (Liu et al., 2014):
Where Nexist is the number of frames which formed targeted hydrogen bond, and Ntotal is the total number of frames. The occupations are varied from 1 to 100%, and a higher percentage represents a more stable-existed hydrogen bond.
Dynamic Cross Correlation Matrix
The dynamic cross-correlation matrix (DCCM) analysis of the Cα atoms during the last 50 ns of the first parallel trajectory was performed to explore the correlated motion between residues of DC3-AR complex. The cross correlation matrix Cij, which reflects the displacements of the Cα atoms relative to average positions, was determined by following equation (Lange et al., 2005; Ghosh and Vishveshwara, 2007):
Where ΔRi and ΔRj represent displacements of atom i and j, respectively. The value of Cij fluctuated from −1 to 1, the positive value indicates a correlated motion between the residue i and residue j, while negative values indicates an anti-correlated motion.
Clustering analysis was conducted by using the MMTSB toolset in Amber 12 to determine the representative structure of DC3-AR complex during the last 50 ns MD simulations. Firstly, the similar conformations of DC3-AR complex generated from the trajectory were classified into one cluster, and the most populated cluster has maximum number of conformations. Centroids of the generated clusters were then calculated and generated. Subsequently, RMSD of each structure was calculated with respect to specific centroid. Ultimately, the structure with lowest RMSD to cluster centroid from the most populated cluster was defined as the representative structure of DC3-AR. After that, the representative structure was adopted to generate the hydrophobic and electrostatic interaction surface of DC3-AR complex by using UCSF Chimera package (Pettersen et al., 2004; http://www.cgl.ucsf.edu/chimera).
Dynamical Correlation Network
The cross-correlation matrix Cij was also employed to build the dynamical correlation network to intuitively exhibit the correlated motion between residues in different protein domains. The Cα atom of each residue was defined as a “node,” and “edge” is the connection of each pair of nodes if the residue pairs interact with each other (Liu et al., 2014). The edges were computed as the following equation:
Here, each edge has a specific contribution to the movement of complex, the motion of residue i can be used to predict the motion direction of residue j. If |i − j|<=10, cross-correlation between i and j are ignored to remove the correlations due to special closeness. Besides, −0.3 ≤ Cij ≤ 0.3 were also deleted to make network plot more concise. Network View plugin in visual Molecular Dynamics 1.9.2 (VMD) (Humphrey et al., 1996; Hsin et al., 2008) was used to visualize the interaction network.
Results and Discussion
Homology Modeling of Cyclopeptide DC3
Through templates searching by SWISS-MODEL, 28 qualified templates for DC3 sequence were found. According to the selection criteria and the basic structure characters of cyclopeptide (Craik et al., 1999; Sze et al., 2009) 2kux.1.A (Plan et al., 2010) was selected as the final template to construct the three-dimensional structure of DC3. As shown in Figure 1A, the amino acid sequence identity between template and target is 56.67%, sequence similarity is 0.52, Global Model Quality Estimation (GMQE) value is 0.94, and the target sequence is all covered. The sequence alignment and structure comparison of target and template were shown in Figure 1B, from where highly similarity can be easily observed. Besides, the Z-score information and predicted local similarity of each residue to target were shown in Figures 1C,D, respectively. All these information reflects the reliability of the constructed model. As the AR residues were numbered from 671 to 919, here, the number of DC3-residues was defined from 641 to 670 for convenience.
Figure 1. Homology modeling results. (A) Information about template, sequence identity, sequence similarity, coverage and GMQE, QMEAN values. (B) Sequence and structure alignment of target and template. The constructed DC3 structure is shown in marine ribbon and the template structure in gray ribbon. G641 and D670 are shown in sticks to exhibit the cyclization state of DC3. (C) Z-score information of the constructed model. (D) The predicted local similarity of each residue to target.
Binding Site Exploration
Molecular Docking Analysis
Molecular docking was then performed to explore the possible binding site and binding mechanism of DC3-AR complex. As a result, 2,000 poses of 60 clusters were generated by ZDOCK module. The pose amount of each cluster, ZRank Score, ZDock Score, and the rationality of binding mode were combined to assess the poses. The ZRank Score represents the extent of energy contribution to the system when a ligand binds to a receptor. The ZDock Score is calculated based on the shape matching degree of receptor and ligand, and a higher score represents a better pose. Furthermore, the electrostatic interaction, Van der Waals' force, and desolvation energy were also taken into consideration. A lower score represents a better pose. Based on these criteria, top four possible binding sites of DC3 to AR were selected. As shown in Figures 2A,B, the top four clusters of DC3-ARs located in Site 1, Site 2, Site 4, and Site 3, respectively. The score of binding to Site 2 was higher than other possible binding sites although the pose number of this cluster is relatively small. However, it should be noted that scoring functions do not always yield the best predictions of binding affinity (Ramírez and Caballero, 2016). To further confirm the binding affinity prediction, Molecular dynamics (MD) simulations were subsequently performed to obtain more conformational sampling of these four systems. Four poses (Table 1, Figures 2C–F) of DC3-AR complex were determined as the initial structures, named as system 1,system 2, system 3, and system 4, to perform MD simulations, which were selected from the top four possible binding sites with good ZRank Score, fine ZDock Score, and rational conformations.
Figure 2. The top four possible binding modes of DC3-AR complex. (A,B) AR is shown in yellow ribbon, and each dot represent a DC3 conformation. The color of dot from red to blue indicates the ZRank Score from high to low. (C–F) Four initial structures of DC3-AR complex selected form the top four possible binding sites for MD simulations. DC3 is shown in marine and AR is shown in violet.
Table 1. The selected top four possible binding sites of DC3-AR complex and the corresponding poses scoring.
Root Mean Square Deviation
One hundred and fifty nano seconds MD simulations were calculated on the four DC3-AR complexes systems acquired from molecular docking were then performed respectively. To obtain reliable and repeatable results, three parallel MD simulations processes were executed on each system. Then root man square deviation (RMSD) values of DC3-AR complexes backbone atoms were calculated relative to the initial structures to monitor the stability and overall convergence of each system during the simulation process. As shown in Figure 3, all systems experienced various degrees of fluctuations at first, but gradually tended to converge. It can be seen that, the first and third trajectories of system 4 and all parallel trajectories of other three systems reached equilibrium in the last 50 ns, which were qualified for subsequent analyses of the dynamic behavior. However, the second parallel trajectory of system 4 experienced great structure changes at about 50 ns, which suggested the relatively poor stability of system 4. Considering the abnormality, this trajectory was eliminated in the succeeding binding free energy analysis which was carried out by averaging the values of parallel trajectories.
Figure 3. Root mean square deviation (RMSD) of backbone atoms of four DC3-AR complexes systems from three parallel trajectories using the initial structure as reference. (A) RMSD of system 1. (B) RMSD of system 2. (C) RMSD of system 3. (D) RMSD of system 4.
Root Mean Square Fluctuation
The root mean square fluctuation (RMSF) reveals the fluctuation of certain residues during simulation process around its average position, which is also a tool to assess the dynamics stability of system. Here, RMSF values of Cα atoms in the last 50 ns were calculated by employing the first parallel trajectory of each system. RMSF of DC3 residues were shown in Figure 4A to explore the stability of the DC3. It can be seen that, different from the great fluctuation in other systems, DC3 residues in system 2 experienced minor motions. It demonstrated that DC3 in system 2 showed obvious superiority in the stability. RMSF of AR residues in these four systems were compared with apo-AR system (made up by AR only) to determine whether the binding of DC3 affects the stability of AR. As shown in Figure 4B, the overall RMSF of system 2 is lower than apo-AR system, especially the residues 840-870 (corresponding to H9, loop 843-849, H10), 880-905 (corresponding to H10, H11, loop888-893, H12). Moreover, only in apo-AR system and system 2, the RMSF values of all residues were under 10Å. Whereas, residues in other systems showed apparently larger conformational changes comparing to apo-AR system. These results demonstrated that the combination of DC3 to AR in site 2 (corresponding to Helix 11, loop 888-893, Helix 12) could stabilize androgen receptor. However, DC3 combination in other sites could visibly reduce AR stability. These RMSF analyses indicated that site 2 is the most possible site of DC3 binding to AR.
Figure 4. Root mean square fluctuation (RMSF) values of Cα atoms calculated by the first parallel trajectory of studied systems in the last 50 ns MD simulation. (A) RMSF values of DC3 residues in four systems. (B) RMSF values of DC3 residues in system2. (C) RMSF values of AR residues in four systems and apo-AR system. (D) RMSF values of AR residues in system 2 and apo-AR system.
Interaction Energetic Features
In order to explore the interaction energetic features of DC3-AR complexes, MM-GBSA method was employed to calculate the binding free energies of each system. The average binding free energies and detailed energetic contribution components of the last 50 ns of parallel trajectories were calculated and shown in Table 2. It can be seen that the free energy of system 2 (−40.94 kcal/mol) is apparently lower than system 1 (−33.49 kcal/mol), system 3 (−18.99 kcal/mol), and system 4 (−19.73 kcal/mol). It demonstrated that DC3 showed a higher binding affinity to AR in system 2 comparing to other systems, which indicated that DC3 has a great tendency to bind to AR in site 2 and system 2 is more likely to remain stable. This result conforms to the conclusion obtained from the previous RMSF analysis. Moreover, details of the dominant components driving DC3 to bind to AR can be acquired by dissecting the binding free energy into contributing components. Here, the electrostatic interaction (ΔEele) in system 2 (−232.64 kcal/mol) can be found to make a great contribution to the low binding free energy of the whole system, which reflects that significant electrostatic interactions may exist between DC3-AR complex and contribute greatly to the system stability.
Table 2. Binding free energy and the detailed energetic contribution components of four systems of DC3-AR complex averaged by the last 50 ns of parallel trajectories (kcal/mol).
Dynamic Cross-Correlation Matrices Analysis
The dynamic cross-correlation matrices (DCCM) analysis was further analyzed (Figure 5) to investigate the correlated conformational motions of DC3-AR complexes. Here, highly positive regions (colored by red and yellow) are associated with strong correlated motions (residue pairs move in the same direction), while negative regions (colored by blue) are linked with strong anti-correlation movements (residue pairs move in the opposite direction). Inspecting the DC3 domains of the four systems, it can be observed that relatively stronger correlations exist between DC3 residues in system 2. Moreover, comparing to other systems, obviously more correlated and anti-correlated motions between DC3 residues and AR residues can be found in system 2. These obvious differences indicated that there were more and stronger cross-correlation motions between residues in system 2, demonstrating more intense interaction and better stability of this DC3-AR complex.
Figure 5. Dynamic cross-correlation matrices (DCCM) of the Cα atoms around their mean positions during the last 50 ns trajectory from the first parallel trajectory. The degrees of correlation and anti-correlation correspond to the color bar. (A) System 1. (B) System 2. (C) System 3. (D) System 4.
Hydrogen Bonds Analysis
The stronger cross-correlation between DC3 and AR residues found in system 2 might also due to the formation of hydrogen bonds during MD simulations. Hydrogen bonds, as critical indicators of nonbonding interactions, play vital roles in the protein-ligand recognition process (Ramírez and Caballero, 2016). During this MD simulation, the number of hydrogen bonds formed between DC3 and AR vs. simulation time was calculated and plotted in Figure 6. Though we set 0.35 nm as the hydrogen bond criteria in this study, the distance we calculated for hydrogen bond were almost all around 3 Angstrom, long-distance hydrogen bond do not exist. As shown in this figure, hydrogen bond interaction patterns formed in system 2 remained constant during the entire simulation time. While in other systems, hydrogen bonds were unstable and most of them disappeared in about 90 ns. Even in the last 50 ns of simulation, the amount of hydrogen bonds still fluctuated a lot. This result reflects the obvious stability of system 2, which further proves that DC3 tends to bind to AR at site 2 as previous RMSF and binding free energy analyses demonstrated. Besides, the intermolecular hydrogen bonds formed in the last 50 ns MD simulations with occupation more than 10% were listed in Table 3. It can be clearly observed that much more hydrogen bonds are stably formed in system 2. On one hand, there are 13 hydrogen bonds occupied more than 10% in system 2, while only 4 hydrogen bonds in system 1, one hydrogen bond in system 3, and even no one in system 4. On the other hand, the highest occupation of hydrogen bonds are 23.51 and 12.55% in system 1 and system 3, respectively, while in system 2 hydrogen bonds formed between K669-H885, D890-R666, S900-T647, and D890-E643 occupied 85.23, 74.54, 72.49, and 64.38%, respectively.
Figure 6. Hydrogen bonds and Distance analyses results. (A) Number of hydrogen bonds formed between DC3 and AR vs. simulation time calculated from the first parallel MD simulation of studied systems. (B) Distance of key residue pairs vs. simulation time calculated from the first parallel MD simulation of system 2.
Table 3. The hydrogen bonds formed between DC3 and AR with occupation more than 10% in the last 50 ns MD simulations of four systems.
Based on these data, it can be concluded that site 2 (H11, loop888-893, H12) is the most possible site of DC3 binding to AR complex.
Exploration of the Binding Mechanism Between DC3 and AR Site 2
To fully explore the binding modes and interaction mechanisms of DC3 and AR, system 2 was further studied to reveal the complicated binding mechanism of DC3-AR complex.
Root Mean Square Fluctuation Analysis
According to the RMSF values of DC3 residues and AR residues shown in Figures 4C,D, respectively, it is easy to determine the key residues dominant in the binding process of DC3-AR complex. In DC3, residues 645-647 and 666 with low RMSF values experienced minor fluctuation, which indicated these residues were relatively more stable during the simulation process. Comparing the AR residues in system 2 to apo-AR system, it can be observed that the RMSF values of residues 840-870 (correspond to H9, loop 843-849, H10), 790-800 (correspond to H7, loop 797-800) in system 2 were lower than apo-AR. which revealed the definite role of these residues in maintaining the system stability. Moreover, residues of binding site 2 (residues 880-903, correspond to H11, loop 888-893, H12) also exhibited obvious lower fluctuation comparing to apo-AR system. It can be observed that the binding site has become one of the most stable regions in system 2. These results not only indicated the dominant role of these residues, but also suggested that specific interactions must have been formed between residues 880-903 and DC3 residues, which then constrained the mobility of them and made the whole system stable.
The representative conformation of DC3-AR complex during MD simulation was extracted by clustering analysis. The first parallel trajectory of system 2 was grouped into 4 clusters based on the conformational similarity. The most populated cluster contained14068 frames, which accounted for 56.27% of all frames extracted from the last 50 ns MD simulation. Then the conformation with least RMSD value in most populated cluster was defined as the representative structure. Here, the representative structure of DC3-AR complex extracted from MD simulations trajectory and the initial structure acquired from molecular docking were plotted in Figures 7A,B, respectively, to exhibit the interaction of key residues visually. From these two figures, it can be seen that residues pairs K669-H885, R666-D890, S643-D890, and R666-E893 were apparently drawn near to each other during the MD simulations. All these indicated that some significant interaction forces might formed between these residues, which further stabilized the complex as RMSF data verified. Similarly, the closeness of residues L654-I882, V657-I882, L658-I882 could also be easily observed in MD simulation through Figures 7C,D. As residues L654, V657, L658, and I882 are nonpolar amino acid, it suggests specific hydrophobic interactions may formed, which needs to be further validated.
Figure 7. Initial structure (A,C) and representative structure (B,D) of the binding site of DC3-AR complex. DC3 is shown in marine and AR is shown in violet. Key residues are shown in stick, and yellow dashed lines represent the distance of specific residue pairs.
Interaction Surface Exploration and Free Energy Decomposition
To figure out the binding mechanism between the key residues in the binding process of DC3-AR complex, hydrophobic and electrostatic interaction surfaces of representative structure were generated. The hydrophobic interaction surface was plotted in Figure 8A, where dodger blue represents hydrophobic minimum, gray depicts the hydrophobicity of 0, and orange represents the largest hydrophobicity. It can be observed that the binding site of AR indeed exits strong hydrophobic interactions with DC3, which promoted the identification and combination of ligand and receptor to a certain extent. Furthermore, two highly hydrophobic interaction domains (deep orange) formed by V649, L650, L651-V901, and L654, V657, L658-L880, L881, I882, respectively, can be found, which played dominant roles in the development of hydrophobic interaction.
Figure 8. Surface of hydrophobic interaction (A) and electrostatic interaction (B) in the binding site of DC3-AR complex. (A) Dodger blue represents hydrophobic minimum, orange represents hydrophobic maximum, key residues are indicated as orange sticks. (B) Positively charged domain colored by blue, negatively charged domain colored by red, key residues are indicated as blue and red sticks.
The electrostatic interaction surface was shown in Figure 8B. It can be seen that most DC3 residues carry positive charge, which come into being a positively charged surface (blue) in the interface. Whereas, a certain number of negatively charged residues (red) existed in binding site of AR. These residues with opposite charges in the interface attracted each other, and made great contribution to the binding process. To deeply investigate the energetic contribution, especially the electrostatic contribution of key residues, free-energy decomposition was performed based on the last 50 ns MD simulations of system 2. The energy contributions of DC3 and AR residues were shown in Figures 9A,B. The electrostatic contributions of DC3 and AR residues were depicted in Figures 9C,D respectively. From the energy contribution it can be seen that residue R666 made incomparable contribution to free-energy, which demonstrated that interactions existed between R666 and AR residues played essential roles in DC3-AR binding process. Meanwhile residues H885, D890, and S900 dominated the energetic contribution. These results proved that residues of both DC3 and AR in binding site do make contributions in the decrease of free energy. Combined with the hydrogen bond analyzed before, it comes to a conclusion that hydrogen bonds formed between K669-H885, D890-R666, S900-T647, and D890-E643 are especially critical components to the interaction between DC3 and AR. Based on all the results above, we can reach a conclusion that K665, R666 of DC3 and E706, E709, D890, E893, E897 of AR ultimately make great contributions to the binding process.
Figure 9. Free-energy decomposition of DC3-AR complex of the first parallel MD simulation in the last 50 ns MD simulation. (A) Energy contribution of DC3 residues. (B) Energy contribution of AR residues in the binding site. (C) Electrostatic contribution of DC3 residues. (D) Electrostatic contribution of AR residues.
To further validate the interaction formed between key residues and investigate the formation process, the distances of key residue pairs mentioned above vs. simulation time were calculated and plotted in Figure 6B. It can be firstly observed that the distance of residue pairs K669-H885, L654-I882, V657-I882, and L658-I882 experienced an obvious decrease in about 40 ns and maintained stable in the later simulation. Significant conformational changes of binding site could be speculated based on this crucial distance variation. This result revalidated the formation of hydrogen bond between K669-H885, and hydrophobic interaction between L645-I882, V657-I882, L658-I882. In addition, residue pairs E643-D890, T647-S900, R666-D890, and R666-E893 kept highly close (about 5Å) and remain stable throughout the simulation process, and some of them even reached about 2.5Å in the last 50 ns. Besides, from the hydrogen bonds figure shown in Figure 6A, it can also be seen that certain hydrogen bonds formed in about 40 ns. All these results proved the stable existence of hydrogen bonds and electrostatic interaction during the whole MD simulations.
To characterize and intuitively exhibit the underlying dynamical cross-correlations among different parts of DC3-AR complex, the overall cross-correlation networks of system 2 and apo-AR system were constructed. As shown in Figure 10A, masses of correlation and anti-correlation both widely and simultaneously existed in apo-AR system, which made it hard to identify the specific cross-correlation pattern, in other words, the interaction between different protein domains was mixed and disorderly. However, when DC3 binding to AR (Figure 10B), all anti-correlation among different AR regions decreased or disappeared, which indicated the improvements of system stability. Besides, organized anti-correlation developed between DC3 and AR regions (loop 687-690, H3, loop 722-725, loop 727-730, H4, loop 822-825, H9, H11, loop 888-893, H12). It reflected the opposite movement tendency between these regions and DC3, namely, they moved close to each other along with the simulation. Moreover, distinct correlation also formed between DC3 and some AR regions (H9, loop 843-849, H10). Combining previous study that residues in these regions had low RMSF values, it can be concluded that this correlation patterns decreased residues fluctuation and enhanced the system stability.
Figure 10. Cross-correlation networks of apo-AR system (A) and DC3-AR system (B). AR is shown in cyan cartoon and DC3 is shown as CPK model. Cross-correlation of residue pairs are distinguished by colors. Correlation are indicated by yellow (0.3 ≤ Cij < 0.5), orange (0.5 ≤ Cij < 0.7) and red (Cij ≥ 0.7). Anti-correlation are indicated by blue (−0.5 < Cij ≤ −0.3) and black (−0.7 < Cij ≤ −0.5).
In this work, the three-dimensional structure of cyclopeptide DC3 was firstly constructed by homology modeling technology using 2kux.1.A as template. Then molecular docking was carried out to predict possible binding site and preferred orientation of DC3 into AR. Finally, four systems with best docking score from top four clusters were selected to perform 150 ns all-atom molecular dynamics (MD) simulations. The MM/GBSA method and a series of MD trajectory analyses were subsequently conducted. The analyses of RMSF, binding free energy, DCCM and hydrogen bonds indicated that DC3 showed a higher binding affinity to AR in site 2 (corresponding to H10, H11, loop888-893, H12) and this system showed obvious superiority in stability comparing to other systems. Besides, much more intermolecular hydrogen bonds were constantly formed in system 2 with high occupation. Stronger cross-correlation among DC3 residues and stronger anti-correlation between DC3 and AR residues also exited here. These results suggest that DC3 is most likely to bind to AR in site 2 encompassed by H10, H11, loop888-893, and H12. Subsequently, combining further analysis of free-energy decomposition, interaction surface, distance, and cross-correlation network, it can be observed that hydrogen bonds, hydrophobic, and electrostatic interactions play dominant roles in the recognition and combination of DC3-AR complex. For hydrogen bonds, it frequently existed between K669-H885, D890-R666, S900-T647, and D890-E643. Besides, K665, R666 of DC3, and E706, E709, D890, E893, E897 of AR made great contributions to electrostatic interaction values. V649, V650, V651-V801, and L654, V657, L658-L880, L881 play essential parts of hydrophobic interaction. These results elucidated the detail interaction mechanism of DC3-AR complex and the key residues dominated in specific interaction. These findings will significantly facilitate our understanding of action mode of DC3 to AR at the molecular level, and contribute to the future rational cyclopeptide drug design for prostate cancer.
JL: conceived and coordinated the study; HZ and TS: did the homology modeling, molecular docking, molecular dynamics simulations, and they wrote the paper; YY and CF: helped to do molecular simulations. All authors analyzed the results and approved the final version of the manuscript.
Conflict of Interest Statement
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.
This research is supported by the Fundamental Research Funds for the Central Universities (lzujbky-2017-203). We would like to thank Prof. Yuguang Mu and Dr. Liangzhen Zheng from School of Biological Sciences of Nanyang Technological University for helping us to treat the data for the figure of cross-correlation network analysis.
Arnold, K., Bordoli, L., Kopp, J., and Schwede, T. (2006). The SWISS-MODEL workspace: a web-based environment for protein structure homology modelling. Bioinformatics 22, 195–201. doi: 10.1093/bioinformatics/bti770
Biasini, M., Bienert, S., Waterhouse, A., Arnold, K., Studer, G., Schmidt, T., et al. (2014). SWISS-MODEL: modelling protein tertiary and quaternary structure using evolutionary information. Nucleic Acids Res. 42, W252–W258. doi: 10.1093/nar/gku340
Bienert, S., Waterhouse, A., de Beer, T. A., Tauriello, G., Studer, G., Bordoli, L., et al. (2016). The SWISS-MODEL repository-new features and functionality. Nucleic Acids Res. 45, D313–D319. doi: 10.1093/nar/gkw1132
Bordoli, L., Kiefer, F., Arnold, K., Benkert, P., Battey, J., and Schwede, T. (2009). Protein structure homology modeling using SWISS-MODEL workspace. Nat. Protoc. 4, 1–13. doi: 10.1038/nprot.2008.197
Center, M. M., Jemal, A., Lortet-Tieulent, J., Ward, E., Ferlay, J., Brawley, O., et al. (2012). International variation in prostate cancer incidence and mortality rates. Eur. Urol. 61, 1079–1092. doi: 10.1016/j.eururo.2012.02.054
Chen, F., Liu, H., Sun, H., Pan, P., Li, Y., Li, D., et al. (2016). Assessing the performance of the MM/PBSA and MM/GBSA methods. 6. Capability to predict protein-protein binding free energies and re-rank binding poses generated by protein-protein docking. Phys. Chem. Chem. Phys. 18, 22129–22139. doi: 10.1039/C6CP03670H
Craik, D. J., Daly, N. L., Bond, T., and Waine, C. (1999). Plant cyclotides: a uniquefamily of cyclic and knotted proteins that defines the cyclic cystine knotstructural motiff. J. Mol. Biol. 294, 1327–1336. doi: 10.1006/jmbi.1999.3383
Darden, T., York, D., Pedersen, L., Darden, T. A., York, D. M., Pedersen, L. G., et al. (1993). Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. doi: 10.1063/1.464397
Fischer, N. M., van Maaren, P. J., Ditz, J. C., Yildirim, A., and van der Spoel, D. (2015). Properties of organic liquids when simulated with long-range Lennard-Jones interactions. J. Chem. Theory Comput. 11, 2938–2944. doi: 10.1021/acs.jctc.5b00190
Fu, T., Wu, X., Xiu, Z., Wang, J., Yin, L., and Li, G. (2013). Understanding the molecular mechanism of binding modes of aurora a inhibitors by long time scale Gpu dynamics. J. Theor. Comput. Chem. 12:1341003. doi: 10.1142/S0219633613410034
Ghosh, A., and Vishveshwara, S. (2007). A study of communication pathways in methionyl- tRNA synthetase by molecular dynamics simulations and structure network analysis. Proc. Natl. Acad. Sci. U.S.A. 104, 15711–15716. doi: 10.1073/pnas.0704459104
Guex, N., Peitsch, M. C., and Schwede, T. (2009). Automated comparative protein structure modeling with SWISS-MODEL and Swiss-PdbViewer: a historical perspective. Electrophoresis 30 (Suppl. 1), S162–S173. doi: 10.1002/elps.200900140
Hitzenberger, M., Schuster, D., and Hofer, T. S. (2017). The binding mode of the sonic hedgehog inhibitor robotnikinin, a combined docking and QM/MM MD study. Front. Chem. 5:76. doi: 10.3389/fchem.2017.00076
Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C. (2006). Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 65, 712–725. doi: 10.1002/prot.21123
Hou, T., Wang, J., Li, Y., and Wang, W. (2011a). 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, 69–82. doi: 10.1021/ci100275a
Hou, T., Wang, J., Li, Y., and Wang, W. (2011b). Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses generated from docking. J. Comput. Chem. 32, 866–877. doi: 10.1002/jcc.21666
Hu, E., Wang, D., Chen, J., and Tao, X. (2015). Novel cyclotides from Hedyotisdiffusa induce apoptosis and inhibit proliferation and migration of prostate cancer cells. Int. J. Clin. Exp. Med. 8, 4059–4065.
Jitonnom, J., Lomthaisong, K., and Lee, V. S. (2012). Computational design of peptide inhibitor based on modifications of proregion from Plutellaxylostella midgut trypsin. Chem. Biol. Drug Des. 79, 583–593. doi: 10.1111/j.1747-0285.2011.01312.x
Kollman, P. A., Massoval, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., et al. (2000). Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. doi: 10.1021/ar000033j
Lange, O. F., Grubmüller, H., and de Groot, B. L. (2005). Molecular dynamics simulations of protein G challenge NMR-derived correlated backbone motions. Angew. Chem. Int. Ed. Engl. 44, 3394–3399. doi: 10.1002/anie.200462957
Lee, H. Z., Bau, D. T., Kuo, C. L., Tsai, R. Y., Chen, Y. C., and Chang, Y. H. (2011). Clarification of the phenotypic characteristics and anti-tumor activity of Hedyotisdiffusa. Am. J. Chin. Med. 39, 201–213. doi: 10.1142/S0192415X11008750
Lin, C. C., Kuo, C. L., Lee, M. H., Hsu, S. C., Huang, A. C., Tang, N. Y., et al. (2011). Extract of Hedyotisdiffusa willd influences murine leukemia WEHI-3 Cells in vivo as well as promoting T- and B-Cell proliferation in leukemic mice. In Vivo 25, 633–640.
Lin, J., Chen, Y., Wei, L., Chen, X., Xu, W., Hong, Z., et al. (2010). Hedyotis Diffusa willd extract induces apoptosis via activation of the mitochondrion-dependent pathway in human colon carcinoma cells. Int. J. Oncol. 37, 1331–1338.
Liu, H., An, X., Li, S., Wang, Y., Li, J., and Liu, H. (2015). Interaction mechanism exploration of R-bicalutamide/S-1 with WT/W741L AR using molecular dynamics simulations. Mol. Biosyst. 11, 3347–3354. doi: 10.1039/C5MB00499C
Liu, M., Wang, L., Sun, X., and Zhao, X. (2014). Investigating the impact of Asp181 point mutations on interactions between PTP1B and phosphotyrosine substrate. Sci. Rep. 4:5095. doi: 10.1038/srep05095
Liu, Z., Liu, M., Liu, M., and Li, J. (2010). Methylanthraquinone from Hedyotisdiffusa willd induces Ca(2+)-mediated apoptosis in human breast cancer cells. Toxicol. In Vitro 24, 142–147. doi: 10.1016/j.tiv.2009.08.002
Mangelsdorf, D. J., Thummel, C., Beato, M., Herrlich, P., Schütz, G., Umesono, K., et al. (1995). The nuclear receptor superfamily: the second decade. Cell 83, 835–839. doi: 10.1016/0092-8674(95)90199-X
Meng, X. Y., Zhang, H. X., Mezei, M., and Cui, M. (2011). Molecular docking: a powerful approach for structure-based drug discovery. Curr. Comput. Aided Drug Des. 7, 146–157. doi: 10.2174/157340911795677602
Onufriev, A., Bashford, D., and Case, D. A. (2004). Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 55, 383–394. doi: 10.1002/prot.20033
Pettersen, E. F., Goddard, T. D., Huang, C. C., Couch, G. S., Greenblatt, D. M., Meng, E. C., et al. (2004). UCSF Chimera–a visualization system for exploratory research and analysis. J. Comput. Chem. 25, 1605–1612. doi: 10.1002/jcc.20084
Plan, M. R., Rosengren, K. J., Sando, L., Daly, N. L., and Craik, D. J. (2010). Structural and biochemical characteristics of the cyclotidekalata B5 from Oldenlandiaaffinis. Biopolymers 94, 647–658. doi: 10.1002/bip.21409
Punkvang, A., Kamsri, P., Saparpakorn, P., Hannongbua, S., Wolschann, P., Irle, S., et al. (2015). Key structures and interactions for binding of mycobacterium tuberculosis protein kinase b inhibitors from molecular dynamics simulation. Chem. Biol. Drug Des. 86, 91–101. doi: 10.1111/cbdd.12465
Ramírez, D., and Caballero, J. (2016). Is it reliable to use common molecular docking methods for comparing the binding affinities of enantiomer pairs for their protein target? Int. J. Mol. Sci. 17:E525. doi: 10.3390/ijms17040525
Ryckaert, J. P., Ciccotti, G., and Berendsen, H. J. C. (1977). Numerical integration of a system with constraints: of the cartesian equations of motion molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. doi: 10.1016/0021-9991(77)90098-5
Schröder, F. H. (2008). Progress in understanding androgen-independent prostate cancer (AIPC): a review of potential endocrine-mediated mechanisms. Eur. Urol. 53, 1129–1137. doi: 10.1016/j.eururo.2008.01.049
Sun, H., Li, Y., Shen, M., Tian, S., Xu, L., Pan, P., et al. (2014a). Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys. 16, 22035–22045. doi: 10.1039/C4CP03179B
Sun, H., Li, Y., Tian, S., Xu, L., and Hou, T. (2014b). Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys. Chem. Chem. Phys. 16, 16719–16729. doi: 10.1039/C4CP01388C
Sze, S. K., Wang, W., Meng, W., Yuan, R., Guo, T., Zhu, Y., et al. (2009). Elucidating the structure of cyclotides by partial acid hydrolysis and LC-MS/MS analysis. Anal. Chem. 81, 1079–1088. doi: 10.1021/ac802175r
Topham, C. M., Thomas, P., Overington, J. P., Johnson, M. S., Eisenmenger, F., and Blundell, T. L. (1990). An assessment of COMPOSER: a rule-based approach to modelling protein structure. Biochem. Soc. Symp. 57, 1–9.
Veldscholte, J., Berrevoets, C. A., Brinkmann, A. O., Grootegoed, J. A., and Mulder, E. (1992). Anti-androgens and the mutated androgen receptor of lncap cells: differential effects on binding affinity, heat-shock protein interaction, and transcription activation. Biochemistry 31, 2393–2399. doi: 10.1021/bi00123a026
Wang, Q., Ning, L., Niu, Y., Liu, H., and Yao, X. (2015). Molecular mechanism of the inhibition and remodeling of human islet amyloid polypeptide (hIAPP(1–37)) oligomer by resveratrol from molecular dynamics simulation. J. Phys. Chem. B 119, 15–24. doi: 10.1021/jp507529f
Wang, Z., Chen, G., Chen, L., Liu, X., Fu, W., Zhang, Y., et al. (2015). Insights into the binding mode of curcumin to MD-2: studies from molecular docking, molecular dynamics simulations and experimental assessments. Mol. Biosyst. 11, 1933–1938. doi: 10.1039/C5MB00085H
Xu, L., Sun, H., Li, Y., Wang, J., and Hou, T. (2013). Assessing the performance of mm/pbsa and mm/gbsa methods. 3. the impact of force fields and ligand charge models. J. Phys. Chem. B 117, 8408–8421. doi: 10.1021/jp404160y
Yamada, A., Fujii, S., Mori, S., and Kagechika, H. (2013). Design and synthesis of 4-(4-Benzoylaminophenoxy)phenol derivatives as androgen receptor antagonists. ACS Med. Chem. Lett. 4, 937–941. doi: 10.1021/ml4001744
Yamamoto, S., Matsunaga, N., Hitaka, T., Yamada, M., Hara, T., Miyazaki, J., et al. (2012). Design, synthesis, and biological evaluation of 4-phenylpyrrole derivatives as novel androgen receptor antagonists. Bioorg. Med. Chem. 20, 422–434. doi: 10.1016/j.bmc.2011.10.067
Yamaoka, M., Hara, T., and Kusaka, M. (2010). Overcoming persistent dependency on androgen signaling after progression to castration-resistant prostate cancer. Clin. Cancer Res. 16, 4319–4324. doi: 10.1158/1078-0432.CCR-10-0255
Zhou, J., Liu, B., Geng, G., and Wu, J. H. (2010). Study of the impact of the T877A mutation on ligand-induced helix-12 positioning of the androgen receptor resulted in design and synthesis of novel antiandrogens. Proteins 78, 623–637. doi: 10.1002/prot.22592
Keywords: Cyclopeptide DC3, androgen receptor, protein drug interaction, homology modeling, molecular docking, molecular dynamics simulations
Citation: Zhang H, Song T, Yang Y, Fu C and Li J (2018) Exploring the Interaction Mechanism Between Cyclopeptide DC3 and Androgen Receptor Using Molecular Dynamics Simulations and Free Energy Calculations. Front. Chem. 6:119. doi: 10.3389/fchem.2018.00119
Received: 29 January 2018; Accepted: 30 March 2018;
Published: 19 April 2018.
Edited by:Sam P. De Visser, University of Manchester, United Kingdom
Reviewed by:Gerardo Andres Cisneros, University of North Texas, United States
Jitrayut Jitonnom, University of Phayao, Thailand
Copyright © 2018 Zhang, Song, Yang, Fu and Li. 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 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: Jiazhong Li, firstname.lastname@example.org
†These authors have contributed equally to this work.