BRIEF RESEARCH REPORT article
Pathways and Mechanism of Caffeine Binding to Human Adenosine A2A Receptor
- Center for Computational Biology and Department of Molecular Biosciences, University of Kansas, Lawrence, KS, United States
Caffeine (CFF) is a common antagonist to the four subtypes of adenosine G-protein-coupled receptors (GPCRs), which are critical drug targets for treating heart failure, cancer, and neurological diseases. However, the pathways and mechanism of CFF binding to the target receptors remain unclear. In this study, we have performed all-atom-enhanced sampling simulations using a robust Gaussian-accelerated molecular dynamics (GaMD) method to elucidate the binding mechanism of CFF to human adenosine A2A receptor (A2AAR). Multiple 500–1,000 ns GaMD simulations captured both binding and dissociation of CFF in the A2AAR. The GaMD-predicted binding poses of CFF were highly consistent with the x-ray crystal conformations with a characteristic hydrogen bond formed between CFF and residue N6.55 in the receptor. In addition, a low-energy intermediate binding conformation was revealed for CFF at the receptor extracellular mouth between ECL2 and TM1. While the ligand-binding pathways of the A2AAR were found similar to those of other class A GPCRs identified from previous studies, the ECL2 with high sequence divergence serves as an attractive target site for designing allosteric modulators as selective drugs of the A2AAR.
Adenosine receptors (ARs) are a subfamily of G-protein-coupled receptors (GPCRs) with adenosine as the endogenous ligands (Fredholm et al., 1997). They belong to class A GPCRs and consist of four known subtypes: A1AR, A2AAR, A2BAR, and A3AR (Jacobson and Gao, 2006). Despite their broad distribution in human tissues and functional differences, ARs share common antagonists of caffeine (CFF) and theophylline, both of which antagonize the receptors upon binding. The sequence alignment by MultiSeq in VMD (Humphrey et al., 1996) showed that the seven transmembrane (TM) helix bundles of the A1AAR shares high similarity with A2AAR by 71%, A2BAR by 70%, and A3AAR by 77%. The sequence similarity is significantly reduced in the three extracellular loops (ECLs), being 43% for A2AAR, 45% for A2BAR, and 35% for A3AR when compared with A1AR. The CFF antagonist binds to all four subtypes of ARs, but with different binding affinities (Porkka-Heiskanen et al., 1997). Understanding the binding mechanism of CFF is expected to facilitate drug design targeting the functionally important ARs.
The human A2AAR is one of the best structurally characterized GPCRs at the atomic level, with more than 30 x-ray and cryo-EM structures published to date (Carpenter and Lebon, 2017). Structures of its three distinct conformational states have been reported, including the inactive conformation bound to an antagonist or inverse agonist, the intermediate conformation bound to an agonist, and the active conformation bound to both an agonist and engineered G protein. The orthosteric binding pocket of the A2AAR is defined by the residues interacting with both the endogenous adenosine agonist and antagonist such as CFF (Carpenter and Lebon, 2017). The residues include M5.38, M7.35, I7.39, and F45.52ECL2 (Doré et al., 2011; Lebon et al., 2011; Cheng et al., 2017). The GPCR residues are numbered according to the Ballesteros–Weinstein scheme (Ballesteros and Weinstein, 1995). Receptor residue N6.55 can form hydrogen bonds with the N7 atom of adenosine (Lebon et al., 2011) and the O11 or O13 atom of CFF, which results in two distinct binding orientations referred to as CFF A and B (Cheng et al., 2017). Other interacting residues include V3.32, L3.33, T3.36, W6.48, L6.51, S7.42, and H7.43 (Lebon et al., 2012; Cheng et al., 2017).
Molecular dynamics (MD) simulations have been previously carried out to characterize binding of the CFF antagonist to the human A2AAR. Cao et al. (2015) performed 800 ns MD simulations to elucidate the effect of membrane composition on the CFF-bound A2AAR. They discovered that the seven TM helix folds were maintained across the systems over the course of their simulations. CFF was flexible and exhibited multiple binding poses in the receptor orthosteric binding pocket. The four most populated binding poses of CFF were extracted with the interacting residues, including A2.61, I2.64, S2.65, V3.32, L3.33, T3.36, F45.52, E169ECL2, M5.38, N5.42, L6.51, H6.52, N6.55, H264ECL3, M7.35, I7.39, and H7.43. In particular, CFF forms a hydrogen bond with receptor residue N6.55 and water-bridge contact with residue H7.43 (Cao et al., 2015). Guo et al. (2016) performed 10 temperature-accelerated MD (TAMD) simulations starting from the 4EIY PDB structure to investigate the dissociation pathway of the ZM241385 antagonist from the A2AAR. The method specifically accelerated the center of mass of the ligand, and thus the A2AAR was almost rigid. They found 16 residues that could potentially interact with ZM241385 during the ligand dissociation process, including G1TM1, I2.63, S2.64, T2.65, Q148ECL2, G152ECL2, K153ECL2, S156ECL2, Q157ECL2, E169ECL2, T6.58, H7.29, A7.30, P7.31, L7.32, and Y7.36. Specifically, the residues E169ECL2, T6.58, and H7.29 along with the structural water of 4EIY formed a hydrogen bond network interacting with the ligand ZM241385 (Guo et al., 2016). Caliman et al. (2018) applied the FTMap fragment-based mapping algorithm on the four distinct conformers obtained from MD simulations of two ligand free receptor conformations of the A2AAR (PDBs: 3QAK and 3EML). They uncovered five non-orthosteric binding sites that were located in the intracellular region of the TM helices TM3/TM4, the G-protein-binding site in the intracellular region between TM2/TM3/TM6/TM7, the lipid interface of TM5/TM6, the intracellular region of TM1/TM7, and the extracellular region of TM3/TM4 of the A2AAR. Their analysis also revealed residues in the orthosteric binding site, including I2.64, V3.32, L3.33, T3.36, Q3.37, I3.40, L45.51ECL2, F45.52ECL2, E169ECL2, M5.38, N5.42, W6.48, L6.51, H6.52, N6.55, T6.58, H264ECL3, L7.32, M7.35, Y7.36, I7.39, S7.42, and H7.43 (Caliman et al., 2018).
Gaussian-accelerated MD (GaMD) is a computational method that allows for simultaneous unconstrained enhanced sampling and free energy calculations of large biomolecules (Miao et al., 2015). By adding a harmonic boost potential, GaMD smooths the potential energy surface of biomolecules to reduce the system energy barriers (Miao et al., 2015). The harmonic boost potential mostly exhibits a Gaussian distribution. Cumulant expansion to the second order (“Gaussian approximation”) can thus be applied to achieve proper energetic reweighting. GaMD resolves the energetic noise problem encountered in the previous accelerated MD (aMD) method (Hamelberg et al., 2004; Shen and Hamelberg, 2008), thereby allowing us to recover the original free energy profiles of biomolecules (Miao et al., 2015). Even though it is exceedingly difficult to obtain convergent free energy profiles for large biomolecular systems, “semiquantitative” low-energy conformational states of biomolecules can be identified from the GaMD-reweighted free energy profiles. GaMD does not require carefully predefined collective variables and as such it is advantageous to study complex biological processes. GaMD has been demonstrated on enhanced sampling and free energy calculations of ligand binding (Miao et al., 2015; Pang et al., 2017), protein folding (Miao et al., 2015; Pang et al., 2017), GPCR activation (Miao and McCammon, 2016), and protein–membrane (Bhattarai et al., 2020), protein–protein (Miao and McCammon, 2018; Wang and Miao, 2019), and protein–nucleic acid (Ricci et al., 2019; East et al., 2020) interactions. Of relevance to studies of GPCRs, GaMD simulations have successfully revealed the mechanisms of GPCR activation, ligand binding, and GPCR–G-protein interactions, which were consistent with experimental data and/or long timescale conventional MD (cMD) simulations (Miao and McCammon, 2016, Miao et al., 2018; Pawnikar and Miao, 2020).
In this study, we have performed all-atom GaMD simulations to determine the pathways and mechanism of CFF binding to the human A2AAR. The GaMD simulations have captured both binding and dissociation of CFF in the A2AAR. The simulation-predicted binding poses were consistent with x-ray crystal conformations of CFF in the 5MZP PDB structure (Cheng et al., 2017). An important intermediate binding site of CFF was also revealed from the GaMD simulations. The simulation findings could provide a molecule basis for rational computer-aided drug design targeting the A2AAR and other ARs.
The x-ray crystal structure of the human A2AAR in complex with CFF at 2.1 Å resolution (PDB: 5MZP) (Cheng et al., 2017) was used for setting up the simulation system. The structure included 296 out of the total 306 residues of the A2AAR, with 10 missing residues (209–218). The T4-lysozyme, lipid molecules, CFF, water, and heteroatom molecules were removed. A total of 10 CFF ligand molecules were placed randomly at a distance >15 Å from the extracellular surface of the A2AAR (Figure 1A). The simulation system was then prepared using the CHARMM-GUI webserver with the membrane input generator (Wang et al., 2006; Jo et al., 2007, 2008; Wu et al., 2014; Lee et al., 2016, 2020). The system dimension was 81.18 × 81.18 × 115.11 Å. It included 161 POPC lipid molecules, with 81 molecules on the upper leaflet and 80 molecules on the lower leaflet, and 14,627 water molecules. All chain termini were capped with neutral patches (acetyl and methylamide). The system was solvated in 0.15 M NaCl solution at temperature 310 K. The AMBER FF19SB (Tian et al., 2019) parameter set was used for the receptor, LIPID17 (Gould et al., in preparation) for the POPC lipids, TIP3P (Jorgensen et al., 1983) for water, and GAFF2 (Wang et al., 2004; He et al., 2020) for CFF. The output files from CHARMM-GUI were used to perform GaMD simulations with AMBER 20 (Case et al., 2020).
Figure 1. Gaussian-accelerated molecular dynamics (GaMD) simulations successfully captured both binding and dissociation of caffeine (CFF) in the A2AAR. (A) Computational model used for simulations of the A2AAR (blue ribbons) with 10 CFF molecules (orange spheres) placed far away in the solvent. The receptor was inserted in a POPC lipid bilayer (cyan sticks) and solvated in an aqueous solution (cyan) of 0.15 M NaCl. (B) X-ray structure of CFF-bound A2AAR (PDB: 5MZP). A hydrogen bond is formed between either O11 or O13 atom of CFF with the ND2 atom of the receptor residue N6.55 in two X-ray conformations of the ligand (CFF-A and CFF-B), in which the distance between the N1 atom that connects atoms O11 and O13 in CFF and the ND2 atom of residue N6.55 stays at 5.1 Å. The seven transmembrane (TM) helices I–VII and three extracellular loops (ECL) 1–3 are labeled in the A2AAR. (C–F) Time courses of the N6.55:ND2–CFF:N1 distance calculated from 63 ns GaMD equilibration and three independent 500–1,000 ns GaMD simulations.
The output files of CHARMM-GUI webserver were used for initial energy minimization, equilibration, and cMD to prepare the system for GaMD simulations. The system was energetically minimized for 5,000 steps using the steepest-descent algorithm and equilibrated with the constant number, volume, and temperature (NVT) ensemble at 310 K using default parameters given by CHARMM-GUI. It was further equilibrated for 375 ps at 310 K with the constant number, pressure, and temperature (NPT) ensemble. The cMD simulation was then performed for 10 ns using the NPT ensemble with constant surface tension at 1 atm pressure and 310 K temperature.
Gaussian-accelerated MD implemented in GPU version of AMBER 20 (Salomon-Ferrer et al., 2013; Miao et al., 2015; Case et al., 2020) was applied to simulate the A2AAR system. The simulations involved an initial short cMD of 3.0 ns to calculate GaMD acceleration parameters and GaMD equilibration of added boost potential for 60 ns. Three independent 500–1,000 ns GaMD production simulations with randomized initial atomic velocities were performed on the A2AAR with 10 unbound CFF molecules. All GaMD simulations were run at the “dual-boost” level by setting the reference energy to the lower bound. One boost potential was applied to the dihedral energetic term and the other to the total potential energetic term. The average and SD of the system potential energies were calculated every 300,000 steps (0.6 ns) for all simulation systems. The upper limit of the boost potential SD, σ0 was set to 6.0 kcal/mol for both the dihedral and the total potential energetic terms. The simulation frames were saved every 1.0 ps for analysis. The GaMD simulations are summarized in Supplementary Table 1.
The NPT ensemble with constant surface tension was used in the short cMD and GaMD simulations. The input files for GaMD equilibration and GaMD simulations have been attached as the Supporting Information. The standard protocol for MD simulations of membrane proteins was followed using notably the system configuration files generated from CHARMM-GUI (Wang et al., 2006; Jo et al., 2007, 2008; Wu et al., 2014; Lee et al., 2016, 2020). Using the MEMBPLUGIN 1.1 plugin of VMD (Guixà-González et al., 2014), we calculated the area per lipid to be 81.58 ± 8.92 Å2 and the membrane thickness to be 70.74 ± 2.09 Å from the GaMD production simulations. The area per lipid was consistent with the initial value of 81.68 Å from CHARMM-GUI. The density of the entire system was calculated to be 1.008 ± 0.001 g/cm3 from the GaMD production simulation outputs, being similar to the value of 1.020 g/cm3 in the cMD simulation. Therefore, the system was expected to behave normally as in other simulation studies.
Simulation analysis was carried out using CPPTRAJ (Roe and Cheatham, 2013) and VMD (Humphrey et al., 1996). The software tools were applied to track the binding and dissociation of CFF from the A2AAR. A hydrogen bond could be formed between O11 or O13 of CFF with atom ND2 in residue N6.55 of the A2AAR, so the distance between atom N1 that connects atoms O11 and O13 in CFF and atom ND2 of residue N6.55 was calculated to monitor ligand binding (Figure 1B). The distance between atom ND2 of receptor residue N6.55 and atom N1 of CFF and the distance of important interactions between receptor residues of TM helices (TM) III, VI, and VII were identified to calculate 2D potential mean force (PMF) free energy profiles using the PyReweighting toolkit (Miao et al., 2014). A bin size of 1 Å was used for the distances. The cutoff was set to 500 frames in one bin for reweighting.
The hierarchical agglomerative clustering algorithm was used to cluster the snapshots of protein conformations with all GaMD production simulations combined. The combined GaMD simulations of CFF binding to the A2AAR were clustered to obtain clusters that corresponded to the low-energy states in the 2D PMF free energy profiles.
Gaussian-Accelerated MD Simulations Captured Both Binding and Dissociation of CFF in the A2AAR
Three independent dual-boost GaMD simulations showed similar averages and SDs of the added boost potentials: 16.21 ± 4.50 kcal/mol for Sim1, 16.20 ± 4.49 kcal/mol for Sim2, and 16.32 ± 4.52 kcal/mol for Sim3, respectively (Supplementary Table 1). Spontaneous binding of CFF to the orthosteric site of the A2AAR was detected at ∼9 ns into the GaMD equilibration (Figure 1C). The first two independent 1,000 ns GaMD simulations (Sim1 and Sim2) captured binding of the CFF in the receptor orthosteric pocket (Figures 1D,E). Remarkably, at ∼400 ns into GaMD Sim2, a second CFF bound to the orthosteric pocket of the A2AAR, while the first CFF remained bound (Figure 1E). The complete dissociation of CFF from the orthosteric pocket of the A2AAR was observed at ∼60 ns in the last independent 500 ns GaMD simulation (Sim3). The traces of CFF binding and dissociation were then analyzed in detail using CPPTRAJ and VMD. The representative CFF poses were selected at distances between receptor residue N6.55 atom ND2 and CFF atom N1 of ∼15, 10, and 5 Å to calculate the interacting residues from the A2AAR in the binding and dissociation pathways using LigPlot (Wallace et al., 1995; Supplementary Figures 1, 2).
Free Energy Profiles of CFF Binding to the A2AAR Receptor
We combined all three GaMD production simulations to calculate reweighted free energy profiles to characterize the binding of CFF to the A2AAR. The distance between the ND2 atom of receptor residue N6.55 and the N1 atom of CFF, ionic lock distance between the CZ atom of residue R3.50 and the CD atom of residue E6.30, and the distance between the CZ atom of residue R3.50 and the OH atom of residue Y7.53 were selected as reaction coordinates to calculate the one-dimensional (1D) (Supplementary Figure 3) and two-dimensional (2D) (Figure 2) PMF free energy profiles. The 1D PMF free energy profiles with variations were obtained by averaging the three GaMD production simulations. Despite the free energy variations, relatively low-energy wells could be identified from 1D PMF profiles of the distances of CFF–residue N6.55, residues R3.50–E6.30, and residues R3.50–Y7.53 (Supplementary Figure 3). The corresponding distances from representative PDB structures (3RFM and 5MZP) were mapped to the 2D free energy profiles for comparison. In the 3RFM PDB structure, the distances between CFF and residue N6.55, residues R3.50 and E6.30, and residues R3.50 and Y7.53 are 5.5, 4.6, and 12.0 Å, respectively. In the 5MZP PDB structure, the distances between CFF and residue N6.55, residues R3.50 and E6.30, and residues R3.50 and Y7.53 are 5.2–5.5, 6.0, and 12.9 Å, respectively.
Figure 2. 2D potential of mean force (PMF) free energy profiles of the A2AAR–caffeine (CFF) interactions. (A) Two-dimensional (2D) PMF of the distance between receptor residue N6.55 atom ND2 and CFF atom N1 and the ionic lock distance between charge centers of receptor residues R3.50 and E6.30. The low-energy states are labeled in the PMF profile, including the unbound (U1 and U2), intermediate (I), and bound (B1 and B2). (B) 2D PMF of the distance between atom CZ in receptor residue R3.50 and the hydroxyl oxygen atom of residue Y7.53 and the ionic lock distance between charge centers of receptor residues R3.50 and E6.30. The low-energy inactive and intermediate conformational states are labeled. The 3RFM and 5MZP PDB structures of inactive A2AAR are mapped to the free energy surface as hexagons and stars.
In the 2D free energy profile of the distances between residue N6.55 and CFF and residues R3.50 and E6.30 (Figure 2A), we identified five low-energy conformational states: unbound (U1, U2), intermediate (I), and bound (B1, B2). In the unbound states (U1, U2), the distance between receptor residue N6.55 and CFF exhibited a broad energy well from ∼35 Å to ∼65 Å, illustrating CFF diffusion in the bulk solvent. The ionic lock distance between residues R3.50 and E6.30 increased from ∼3.5–4.5 Å in U1 to ∼6–7.5 Å in U2. The intermediate state (I) was identified at ∼20–25 Å distance between receptor residue N6.55 and CFF and at ∼7 Å distance between residues R3.50 and E6.30, suggesting that CFF was located at the extracellular mouth of the A2AAR. Both the bound states (B1, B2) were observed at ∼5–10 Å distance between receptor residue N6.55 and CFF. CFF was located in the orthosteric pocket in these states. Similar to the unbound states, the ionic lock distance was ∼3.5–4.5 Å in B1 and ∼6.5–7 Å in B2.
We identified two low-energy conformational states from the free energy profile of the R3.50–Y7.53 and R3.50–E6.30 distances in Figure 2B, labeled as the inactive and intermediate states. In the inactive state, the ionic lock distance between receptor residues R3.50 and E6.30 was ∼3.5–4.5 Å and the distance between receptor residues R3.50 and Y7.53 was ∼10–15 Å. In the intermediate state, the ionic lock distance remained the same, but the distance between receptor residues R3.50 and Y7.53 decreased to ∼5–8 Å.
Binding and Dissociation Pathways of CFF in the A2AAR
In the equilibration trajectory of GaMD simulation, 1 out of the 10 CFF molecules (CFF2) that freely diffused in the solvent bound to the A2AAR through a pathway connecting ECL2, the extracellular mouth between ECL2 and ECL3, and finally the receptor orthosteric site (Figure 3A and Supplementary Figure 1A). At ∼15 Å distance between CFF and the receptor residue N6.55, CFF interacted with the receptor N-terminus and ECL2 (Supplementary Figure 1B). At ∼10 Å distance between CFF and the receptor residue N6.55, CFF was located at the extracellular mouth of the A2AAR between ECL2, ECL3, and TM6, interacting with residues L45.51ECL2, E169ECL2, S263ECL3, and T6.58 (Supplementary Figure 1C). At ∼5 Å distance between CFF and receptor residue N6.55, CFF bound to the receptor orthosteric site. In Sim2 of GaMD simulation trajectory, out of the nine remaining CFF molecules that freely diffused in the solvent, another CFF (CFF9) bound to the orthosteric pocket of the A2AAR, while CFF2 remained bound (Figures 1E, 3B). The binding pathway of CFF9 was mostly similar to that of CFF2, except a slight difference that CFF9 explored a region between ECL2 and TM6 after entry into the receptor (Figure 3B).
Figure 3. Binding and dissociation pathways of caffeine (CFF) in the A2AAR revealed from the Gaussian-accelerated molecular dynamics (GaMD) simulations. (A) Trace of CFF (orange) binding to the A2AAR observed in the GaMD equilibration. Starting from free diffusion in the solvent, CFF bound to the orthosteric site of the A2AAR receptor. (B) Binding of the second CFF (orange) to orthosteric pocket of the A2AAR observed in GaMD Sim2. The first bound CFF is shown in red. (C) Pathway of CFF that dissociated from orthosteric site of the A2AAR to the bulk solvent observed in GaMD Sim3. The A2AAR receptor is shown in blue ribbons and the CFF traces are shown as orange beads. (D) The B1-bound conformational state of CFF was located between ECL2, TM3, TM5, and TM6 with interacting residues F45.52ECL2, V3.32, M5.38, N5.42, W6.48, L6.51, H6.52, and N6.55. (E) The B2-bound conformational state of CFF was located between ECL2, TM3, TM5, and TM6 with interacting residues V3.32, L3.33, F45.52ECL2, M5.38, N5.42, and N6.55. CFF formed a hydrogen bond with the receptor residue N6.55 in both the B1 and B2 states. (F) The intermediate (I) conformational state of CFF that was located between ECL2, N-terminus of TM1, and TM2 with interacting residues P1.28, I1.29, S2.65, K153ECL2, S156ECL2, Q157ECL2, and L45.51ECL2. CFF formed hydrogen bonds with both receptor residues I1.29 and Q157ECL2. The CFF ligand is represented by sticks with carbon atoms colored in orange for simulation-derived low-energy conformations and pink and purple for two x-ray conformations in the 5MZP PDB structure. The receptor-interacting residues are highlighted in green.
In Sim3 of GaMD simulation trajectory, CFF dissociated from the orthosteric site of the A2AAR to the bulk solvent through a pathway connecting the receptor orthosteric pocket and the extracellular mouth between ECL2 and TM7 (Figure 3C and Supplementary Figure 2A). At ∼5 Å distance between CFF and the receptor residue N6.55, CFF bound to the orthosteric site of the A2AAR (Supplementary Figure 2B). At ∼10 Å distance between CFF and the receptor residue N6.55, CFF was located at the extracellular mouth of the A2AAR between ECL2 and TM7, interacting with residues L45.51ECL2, I1.29, I2.64, L7.32, and Y7.36 (Supplementary Figure 2C). At ∼15 Å distance between CFF and the receptor residue N6.55, CFF moved near ECL2–TM1 and interacted with receptor residues L45.51ECL2, A1.26, P1.27, P1.28, and I1.29 (Supplementary Figure 2D). At ∼20 Å distance between CFF and the receptor residue N6.55, CFF is in the intermediate (I) conformational state. While the GaMD simulations were not sufficiently converged with only a few ligand-binding events captured, the binding and dissociation pathways of CFF characterized using the GaMD energetically reweighted structural clusters of the ligand (Supplementary Figure 4) were similar to those as shown in Figure 2.
Low-Energy Binding Poses of CFF in the A2AAR
Next, we combined GaMD simulations of CFF binding to the A2AAR and clustered the simulation snapshots of CFF to obtain representative structural clusters that corresponded to the low-energy states in the 2D PMF free energy profiles (Figure 2A). In the B1-bound state, CFF bound to the orthosteric pocket of the A2AAR and interacted with residues F45.52ECL2, V3.32, M5.38, N5.42, W6.48, L6.51, H6.52, and N6.55. In particular, a hydrogen bond was formed between the ND2 atom of receptor residue N6.55 and O13 atom of CFF at a distance of 2.7 Å (Figure 3D and Supplementary Figure 1D). In the B2-bound state, CFF bound to the orthosteric pocket of the A2AAR in the presence of another CFF molecule in the pocket. The orthosteric pocket was located within the receptor TM bundle between ECL2, TM3, TM5, and TM6. CFF interacted with residues F45.52ECL2, V3.32, L3.33, M5.38, N5.42, and N6.55. In particular, a hydrogen bond was formed between the ND2 atom of receptor residue N6.55 and O13 atom of CFF at a distance of 2.9 Å (Figure 3E).
In the intermediate (I) conformational state (Figure 3F), CFF was located at the extracellular mouth of the A2AAR between ECL2 and TM1, interacting with residues P1.28, I1.29, S2.65, K153ECL2, S156ECL2, Q157ECL2, and L45.51ECL2. In particular, a hydrogen bond was formed between the N atom of receptor residue I1.29 and N9 atom of CFF and another hydrogen bond was formed between the NE2 atom of receptor residue Q157ECL2 and O13 atom of CFF (Figure 3F).
Conformational Changes of the A2AAR During CFF Binding
Two different low-energy conformational states were identified from GaMD simulations of the A2AAR during CFF binding, including the inactive and intermediate states (Figure 2B). The hierarchical agglomerative clustering algorithm was used to cluster snapshots of the A2AAR conformations with all the GaMD production simulations combined. The combined GaMD simulation trajectories were clustered to identify representative low-energy conformational states of the receptor (Figure 4). The ionic lock distance between residues R3.50 and E6.30 changed from 6 Å in the 5MZP PDB structure to 4.3–4.4 Å in the inactive and intermediate conformations (Figure 4A). The distance between the atom CZ of residue R3.50 and the OH atom of residue Y7.53 decreased from 12.4 Å in the inactive conformation (similar to 13.2 Å in the 5MZP PDB structure) to 6.0 Å in the intermediate cluster (Figure 4B). Therefore, the NPxxY motif, a highly conserved motif in the intracellular end of TM7 of class A GPCRs, moved inward during the conformational transition of the A2AAR from the inactive to the intermediate state.
Figure 4. Representative “inactive” (green) and “intermediate” (orange) low-energy conformations of the A2AAR compared with the 5MZP PDB structure (blue). (A) The ionic lock distance between charge centers of residues R3.50 and E6.30 in the 5MZP, inactive, and intermediate conformations are 6.0, 4.4, and 4.3 Å, respectively. (B) The distance between atom CZ in residue R3.50 and hydroxyl oxygen atom of Y7.53 in the 5MZP, inactive, and intermediate conformations are 13.2, 12.4, and 6.0 Å, respectively.
In this study, all-atom GaMD simulations have been applied to elucidate the pathways and mechanism of CFF binding to the human A2AAR. The GaMD simulations have successfully captured both spontaneous binding and dissociation of CFF in the receptor. With GaMD-enhanced sampling, we were able to simulate the complete binding of the CFF antagonist with the final orthosteric pocket deeply buried in the receptor TM domain. However, it is important to note that only two ligand-binding events and one dissociation event were observed in the presented GaMD simulations (Figure 1). Quantitative characterization of the ligand-binding free energy and kinetics would require sampling of significantly more ligand-binding events, which will be investigated in the future using a very recently developed, potentially more efficient Ligand GaMD (LiGaMD) method (Miao et al., 2020) and other applicable algorithms. Nevertheless, energetic reweighting of the GaMD simulations enabled us to identify relatively low-energy conformational states of CFF binding to the A2AAR. Our results were consistent with the experimental data from 3RFM PDB (Doré et al., 2011) and 5MZP PDB (Cheng et al., 2017). In the 3RFM PDB structure, the residues interacting with CFF include F45.52ECL2, M5.38, L6.51, N6.55, M7.35, and I7.39. In the 5MZP PDB structure, the residues interacting with CFF include I2.64, V3.32, F45.52ECL2, M5.38, L6.51, N6.55, and I7.39. The B1- and B2-bound conformations identified from the GaMD free energy profiles were comparable with the experimental 3RFM and 5MZP PDB structures in terms of the ligand interacting residues in the orthosteric pocket and distances between residues R3.50 and E6.30 and residues R3.50 and Y7.53.
We identified a dominant pathway of CFF binding to the A2AAR from the GaMD simulations. CFF approached the A2AAR through interactions with ECL2, extracellular mouth between ECL2, ECL3, and TM7, and finally the receptor orthosteric site located deeply within the receptor TM bundle (Figure 3A and Supplementary Figure 1). A slightly different binding pathway was observed when two CFF molecules bound to the orthosteric pocket of the A2AAR. In this pathway, the second CFF explored a region between ECL3 and TM7 during the binding process (Figure 3B). The dissociation pathway of CFF observed from the GaMD simulation was mostly the reverse of the dominant binding pathway (Figure 2 and Supplementary Figure 4). CFF moved from the receptor orthosteric site to the extracellular mouth between ECL2 and TM7 and then ECL2 and TM1 before dissociating to the bulk solvent (Figure 3C and Supplementary Figure 2). Two low-energy conformational states were identified from the GaMD simulations of the A2AAR during CFF binding, i.e., the inactive and intermediate states. In the inactive state, the distances between residues R3.50 and E6.30 and R3.50 and Y7.53 were 4.4 and 12.4 Å, respectively. In this context, the average distances between residues R3.50 and E6.30 and R3.50 and Y7.53 in 46 experimental structures of the inactive A2AAR (Supplementary Table 2; Jaakola et al., 2008; Doré et al., 2011; Congreve et al., 2012; Hino et al., 2012; Liu et al., 2012; Batyuk et al., 2016; Segala et al., 2016; Cheng et al., 2017; Martin-Garcia et al., 2017, 2019; Melnikov et al., 2017; Sun et al., 2017; Weinert et al., 2017; Broecker et al., 2018; Eddy et al., 2018; Rucktooa et al., 2018; Ishchenko et al., 2019; Shimazu et al., 2019; Borodovsky et al., 2020; Ihara et al., 2020; Jespers et al., 2020; Lee et al., 2020; Nass et al., 2020) were calculated to be 6.5 ± 1.0Å and 12.7 ± 0.4Å, respectively. Therefore, the highly conserved residues R3.50 and E6.30 ionic lock became fully closed in the GaMD simulations of the inactive A2AAR during binding of the CFF antagonist. In comparison, the average distances between residues R3.50 and E6.30 and residues R3.50 and Y7.53 in the nine available structures of active A2AAR (Supplementary Table 2; Lebon et al., 2011, 2015; Xu et al., 2011; Carpenter et al., 2016; Garcia-Nafria et al., 2018; White et al., 2018) were calculated as 11.1 ± 0.4Å and 4.4 ± 0.2Å, respectively. No intermediate structure is currently available for the A2AAR (Pándy-Szekeres et al., 2017). In the GaMD-predicted intermediate conformational state of the A2AAR, the ionic lock distance between residues R3.50 and E6.30 was 4.3 Å, similar to that in the inactive receptor. The distance between residues R3.50 and Y7.53, however, decreased to 6.0 Å, comparable with the average of active A2AAR structures. Therefore, while the ionic lock remained closed, the conserved NPxxY motif in the intracellular end of TM7 was able to move inward in the intermediate state of the A2AAR, being consistent with the pathway and mechanism of GPCR activation revealed from earlier studies (Dror et al., 2011; Miao et al., 2013).
The binding and dissociation of CFF antagonist in our GaMD simulations of the A2AAR involved receptor residues P1.28, I1.29, S2.65, L45.51ECL2, F45.52ECL2, K153ECL2, S156ECL2, Q157ECL2, E169ECL2, A259ECL3, S263ECL3, H264ECL3, N6.55, T6.58, F6.59, P7.31, L7.32, L7.34, M7.35, and Y7.36 (Figure 3 and Supplementary Figures 1, 2). The orthosteric pocket was located within the receptor TM bundle and made of receptor residues F45.52ECL2, V3.32, L3.33, M5.38, N5.42, W6.48, L6.51, H6.52, and N6.55. Notably, CFF formed a hydrogen bond with receptor residue N6.55. The four most populated CFF binding poses in the A2AAR found by Cao et al. consisted of residues A2.61, I2.64, S2.65, V3.32, L3.33, T3.36, F45.52, E169ECL2, M5.38, N5.42, L6.51, H6.52, N6.55, H264ECL3, M7.35, I7.39, and H7.43. Furthermore, CFF formed a hydrogen bond with receptor residue N6.55 and water-bridge contact with residue H7.43 (Cao et al., 2015). The ligand dissociation pathway in the A2AAR discovered by Guo et al. involved 16 receptor residues: G1TM1, I2.63, S2.64, T2.65, Q148ECL2, G152ECL2, K153ECL2, S156ECL2, Q157ECL2, E169ECL2, T6.58, H7.29, A7.30, P7.31, L7.32, and Y7.36 (Guo et al., 2016). The residues in the orthosteric binding site of the A2AAR revealed by Caliman et al. (2018) were I2.64, V3.32, L3.33, T3.36, Q3.37, I3.40, L45.51ECL2, F45.52ECL2, E169ECL2, M5.38, N5.42, W6.48, L6.51, H6.52, N6.55, T6.58, H264ECL3, L7.32, M7.35, Y7.36, I7.39, S7.42, and H7.43. Overall, our results were in good agreement with previous studies of the A2AAR, in terms of the receptor residues involved in the ligand dissociation and binding.
An intermediate ligand-binding site was also revealed from the GaMD simulations of CFF binding and dissociation in the A2AAR. It was located at the extracellular mouth between ECL2 and TM1 of the A2AAR. This region has been identified as an allosteric site of many class A GPCRs (Dror et al., 2013; Kruse et al., 2013; Miao and McCammon, 2016; Miao et al., 2018; Pawnikar and Miao, 2020). Taken together, our simulations suggest that CFF binds to the orthosteric pocket of A2AAR via an intermediate site located at the receptor extracellular mouth. The ECL2 with high sequence divergence could serve as an attractive target site for designing allosteric modulators as selective drugs of the A2AAR and other ARs (Miao et al., 2018).
Data Availability Statement
The authors acknowledge that the data presented in this study must be deposited and made publicly available in an acceptable repository, prior to publication. Frontiers cannot accept a manuscript that does not adhere to our open data policies.
YM designed the research. HD performed the research. HD, SA, and YM analyzed the data. HD and YM wrote the manuscript. All authors contributed to the article and approved the submitted version.
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.
We appreciate the preliminary simulations of ARs by Keeley Collins and Amar Kumar. This work used supercomputing resources with allocation award TG-MCB180049 through the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562, and project M2874 through the National Energy Research Scientific Computing Center (NERSC), which is a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231, and the Research Computing Cluster at the University of Kansas. This work was supported in part by the American Heart Association (Award 17SDG33370094), the National Institutes of Health (R01GM132572), and the startup funding in the College of Liberal Arts and Sciences at the University of Kansas.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.673170/full#supplementary-material
Ballesteros, J. A., and Weinstein, H. (1995). Integrated methods for the construction of three-dimensional models and computational probing of structure-function relations in G protein-coupled receptors. Methods Neurosci. 25, 366–428. doi: 10.1016/s1043-9471(05)80049-7
Batyuk, A., Galli, L., Ishchenko, A., Han, G. W., Gati, C., Popov, P. A., et al. (2016). Native phasing of x-ray free-electron laser data for a G protein-coupled receptor. Sci. Adv. 2:e1600292. doi: 10.1126/sciadv.1600292
Borodovsky, A., Barbon, C. M., Wang, Y., Ye, M., Prickett, L., Chandra, D., et al. (2020). Small molecule AZD4635 inhibitor of A2AR signaling rescues immune cell function including CD103+ dendritic cells enhancing anti-tumor immunity. J. Immunother. Cancer 8:e000417. doi: 10.1136/jitc-2019-000417
Broecker, J., Morizumi, T., Ou, W.-L., Klingel, V., Kuo, A., Kissick, D. J., et al. (2018). High-throughput in situ X-ray screening of and data collection from protein crystals at room temperature and under cryogenic conditions. Nat. Protoc. 13, 260–292. doi: 10.1038/nprot.2017.135
Cao, R. Y., Rossetti, G., Bauer, A., and Carloni, P. (2015). Binding of the antagonist CFF to the human adenosine receptor hA2AR in nearly physiological conditions. PLoS One 10:e0126833. doi: 10.1371/journal.pone.0126833
Cheng, R. K. Y., Segala, E., Robertson, N., Deflorian, F., Doré, A. S., Errey, J. C., et al. (2017). Structures of human A1 and A2A adenosine receptors with xanthines reveal determinants of selectivity. Structure 25, 1275–1285.e4.
Congreve, M., Andrews, S. P., Doré, A. S., Hollenstein, K., Hurrell, E., Langmead, C. J., et al. (2012). Discovery of 1,2,4-Triazine derivatives as adenosine A(2A) antagonists using structure based drug design. J. Med. Chem. 55, 1898–1903. doi: 10.1021/jm201376w
Doré, A. S., Robertson, N., Errey, J. C., Ng, I., Hollenstein, K., Tehan, B., et al. (2011). Structure of the adenosine A2A receptor in complex with ZM241385 and the Xanthines XAC and CFF. Structure 19, 1283–1293. doi: 10.1016/j.str.2011.06.014
Dror, R. O., Green, H. F., Valant, C., Borhani, D. W., Valcourt, J. R., Pan, A. C., et al. (2013). Structural basis for modulation of a G-protein coupled receptor by allosteric drugs. Nature 503, 295–299. doi: 10.1038/nature12595
East, K. W., Newton, J. C., Morzan, U. N., Narkhede, Y. B., Acharya, A., Skeens, E., et al. (2020). Allosteric motions of the CRISPR-Cas9 HNH nuclease probed by NMR and molecular dynamics. J. Am. Chem. Soc. 142, 1348–1358. doi: 10.1021/jacs.9b10521
Eddy, M., Lee, M.-Y., Gao, Z.-G., White, K. L., Didenko, T., Horst, R., et al. (2018). Allosteric coupling of drug binding and intracellular signaling in the A2A adenosine receptor. Cell 172, 68–80.e12.
Fredholm, B., Abbracchio, M. P., Burnstock, G., Dubyak, G. R., Harden, T. K., Jacobson, K. A., et al. (1997). Towards a revised nomenclature for P1 and P2 receptors. Trends Pharmacol. Sci. 18, 79–82. doi: 10.1016/s0165-6147(96)01038-3
Guixà-González, R., Rodriguez-Espigares, I., Ramíez-Anguita, J. M., Carrió-Gaspar, P., Martinez-Seara, H., Giorgino, T., et al. (2014). MEMBPLUGIN: studying membrane complexity in VMD. Bioinformatics 30, 1478–1480. doi: 10.1093/bioinformatics/btu037
Guo, D., Pan, A. C., Dror, R. O., Mocking, T., Liu, R., Heitman, L. H., et al. (2016). Molecular basis of ligand dissociation from the adenosine A2A receptor. Mol. Pharmacol. 89, 485–491. doi: 10.1124/mol.115.102657
Hamelberg, D., Mongan, J., and McCammon, J. (2004). Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules. J. Chem. Phys. 120, 11919–11929. doi: 10.1063/1.1755656
He, X., Man, V. H., Yang, W., Lee, T.-S., and Wang, J. (2020). A fast and high-quality charge model for the next generation general AMBER force field. J. Chem. Phys. 153:114502. doi: 10.1063/5.0019056
Hino, T., Arakawa, T., Iwanari, H., Yurugi-Kobayashi, T., Ikeda-Suno, C., Nakada-Nakura, Y., et al. (2012). G-protein-coupled receptor inactivation by an allosteric inverse-agonist antibody. Nature 482, 237–240. doi: 10.1038/nature10750
Ihara, K., Hato, M., Nakane, T., Yamashita, K., Kimura-Someya, T., Hosaka, T., et al. (2020). Isoprenoid-chained lipid EROCOC 17+4 : a new matrix for membrane protein crystallization and a crystal delivery medium in serial femtosecond crystallography. Sci. Rep. 10:19305.
Ishchenko, A., Stauch, B., Han, G. W., Batyuk, A., Shiriaeva, A., Li, C., et al. (2019). Toward G protein-coupled receptor structure-based drug design using X-ray lasers. IUCrJ 6, 1106–1119. doi: 10.1107/s2052252519013137
Jaakola, V.-P., Griffith, M. T., Hanson, M. A., Cherezov, V., Chien, E. Y. T., Lane, J. R., et al. (2008). The 2.6 angstrom crystal structure of a human A2A adenosine receptor bound to an antagonist. Science 322, 1211–1217. doi: 10.1126/science.1164772
Jespers, W., Verdon, G., Azuaje, J., Majellaro, M., Keranen, H., García-Mera, X., et al. (2020). X-Ray crystallography and free energy calculations reveal the binding mechanism of A 2A adenosine receptor antagonists. Angew. Chem. Int. Ed. Engl. 59, 16536–16543. doi: 10.1002/anie.202003788
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi: 10.1063/1.445869
Lebon, G., Warne, T., Edwards, P. C., Bennett, K., Langmead, C. J., Leslie, A. G. W., et al. (2011). Agonist-bound adenosine A2A receptor structures reveal common features of GPCR activation. Nature 474, 521–525. doi: 10.1038/nature10136
Lee, J., Cheng, X., Swalis, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., et al. (2016). CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. J. Chem. Theory Comput. 12, 405–413. doi: 10.1021/acs.jctc.5b00935
Lee, M.-Y., Geiger, J., Ishchenko, A., Han, G. W., Barty, A., White, T. A., et al. (2020). Harnessing the power of an X-ray laser for serial crystallography of membrane proteins crystallized in lipidic cubic phase. IUCrJ 7, 976–984. doi: 10.1107/s2052252520012701
Liu, W., Chun, E., Thompson, A. A., Chubukov, P., Xu, F., Katritch, V., et al. (2012). Structural basis for allosteric regulation of GPCRs by sodium ions. Science 337, 232–236. doi: 10.1126/science.1219218
Martin-Garcia, J., Conrad, C. E., Nelson, G., Stander, N., Zatsepin, N. A., Zook, J., et al. (2017). Serial millisecond crystallography of membrane and soluble protein microcrystals using synchrotron radiation. IUCrJ 4, 439–454.
Martin-Garcia, J., Zhu, L., Mendez, D., Lee, M.-Y., Chun, E., Li, C., et al. (2019). High-viscosity injector-based pink-beam serial crystallography of microcrystals at a synchrotron radiation source. IUCrJ 6, 412–425. doi: 10.1107/s205225251900263x
Melnikov, I., Polovinkin, V., Kovalev, K., Gushchin, I., Shevtsov, M., Shevchenko, V., et al. (2017). Fast iodide-SAD phasing for high-throughput membrane protein structure determination. Sci. Adv. 3:e1602952. doi: 10.1126/sciadv.1602952
Miao, Y., Bhattarai, A., and Wang, J. (2020). Ligand gaussian accelerated molecular dynamics (LiGaMD): characterization of ligand binding and thermodynamics and kinetics. J. Chem. Theory Comput. 16, 5526–5547.
Miao, Y., Feher, V. A., and McCammon, J. A. (2015). Gaussian accelerated molecular dynamics: unconstrained enhanced sampling and free energy calculation. J. Chem. Theory Comput. 11, 3584–3595. doi: 10.1021/acs.jctc.5b00436
Miao, Y., and McCammon, J. A. (2016). Graded activation and free energy landscapes of a muscarinic G-protein-coupled receptor. Proc. Natl. Acad. Sci. U.S.A. 113, 12162–12167. doi: 10.1073/pnas.1614538113
Miao, Y., and McCammon, J. A. (2018). Mechanism of the G-protein mimetic nanobody binding to a muscarinic G-protein-coupled receptor. Proc. Natl. Acad. Sci. U.S.A. 115, 3036–3041. doi: 10.1073/pnas.1800756115
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, 2677–2689. doi: 10.1021/ct500090q
Pawnikar, S., and Miao, Y. (2020). Pathway and mechanism of drug binding to chemokine receptors revealed by accelerated molecular simulations. Future Med. Chem. 12, 1213–1225. doi: 10.4155/fmc-2020-0044
Porkka-Heiskanen, T., Strecker, R. E., Thakkar, M., Bjorkum, A. A., Greene, R. W., and McCarley, R. W. (1997). Adenosine: a mediator of the sleep-inducing effects of prolonged wakefulness. Science 276, 1265–1268. doi: 10.1126/science.276.5316.1265
Ricci, C. G., Chen, J. S., Miao, Y., Jinek, M., Doudna, J. A., McCammon, J. A., et al. (2019). Deciphering Off-target effects in CRISPR-Cas9 through accelerated molecular dynamics. ACS Cent. Sci. 5, 651–662. doi: 10.1021/acscentsci.9b00020
Rucktooa, P., Cheng, R. K. Y., Segala, E., Geng, T., Errey, J. C., Brown, G. A., et al. (2018). Towards high throughput GPCR crystallography: In Meso soaking of Adenosine A2A receptor crystals. Sci. Rep. 8:41.
Salomon-Ferrer, R., Gotz, A. W., Poole, D., Le Grand, S., and Walker, R. C. (2013). Routined microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent Particle Mesh Ewald. J. Chem. Theory Comput. 9, 3878–3888. doi: 10.1021/ct400314y
Segala, E., Guo, D., Cheng, R. K. Y., Bortolato, A., Deflorian, F., Doré, A. S., et al. (2016). Controlling the dissociation of ligands from the adenosine A2A receptor through modulation of salt bridge strength. J. Med. Chem. 59, 6470–6479. doi: 10.1021/acs.jmedchem.6b00653
Shimazu, Y., Tono, K., Tanaka, T., Yamanaka, Y., Nakane, T., Mori, C., et al. (2019). High-viscosity sample-injection device for serial femtosecond crystallography at atmospheric pressure. J. Appl. Crystallogr. 52, 1280–1288. doi: 10.1107/s1600576719012846
Sun, B., Bachhawat, P., Chu, M. L.-H., Wood, M., Ceska, T., Sands, Z. A., et al. (2017). Crystal structure of the adenosine A2A receptor bound to an antagonist reveals a potential allosteric pocket. Proc. Natl. Acad. Sci. U.S.A. 114, 2066–2071. doi: 10.1073/pnas.1621423114
Tian, C., Kasavajhala, K., Belfon, K. A. A., Raguette, L., Huang, H., Migues, A. N., et al. (2019). ff19SB: amino-acid specific protein backbone parameters trained against quantum mechanics energy surfaces in solution. J. Chem. Theory Comput. 16, 528–552. doi: 10.1021/acs.jctc.9b00591
Wallace, A., Laskowski, R., and Thornton, J. (1995). LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions. Protein Eng. Des. Select. 8, 127–134. doi: 10.1093/protein/8.2.127
Weinert, T., Olieric, N., Cheng, R., Brunle, S., James, D., Ozerov, D., et al. (2017). Serial millisecond crystallography for routine room-temperature structure determination at synchrotrons. Nat. Commun. 8: 542.
White, K., Eddy, M., Gao, Z.-G., Han, G. W., Lian, T., Deary, A., et al. (2018). Structural connection between activation microswitch and allosteric sodium site in GPCR signaling. Structure 26, 259–269.e5.
Keywords: adenosine A2A receptor, caffeine, Gaussian accelerated molecular dynamics, ligand binding, mechanism, pathways
Citation: Do HN, Akhter S and Miao Y (2021) Pathways and Mechanism of Caffeine Binding to Human Adenosine A2A Receptor. Front. Mol. Biosci. 8:673170. doi: 10.3389/fmolb.2021.673170
Received: 26 February 2021; Accepted: 24 March 2021;
Published: 27 April 2021.
Edited by:Yong Wang, University of Copenhagen, Denmark
Reviewed by:Lei Fu, Beijing Normal University, China
Haohao Fu, Nankai University, China
James Joseph McCarty, Western Washington University, United States
Copyright © 2021 Do, Akhter and Miao. 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: Yinglong Miao, email@example.com