Yue Pan1, Renrui Qi1, Minghao Li1, Bingda Wang1, Honglan Huang2, Weiwei Han1. 1. Key Laboratory for Molecular Enzymology and Engineering of Ministry of Education, School of Life Science, Jilin University 2699 Qianjin Street Changchun 130012 China weiweihan@jlu.edu.cn. 2. Department of Pathobiology, College of Basic Medical Sciences, Jilin University Changchun China hhl@jlu.edu.cn.
Adenosine deaminase (ADA; EC 3.5.4.4) is an important enzyme in purine metabolism. ADA is widely distributed in various tissues of the human body and catalyzes the irreversible deamination of adenosine to inosine.[1-5] ADA deficiency is related to various immune deficiency (SCID) diseases with profound depletion of T, B, and natural killer cell lineages, such as rheumatoid arthritis, an autoimmune disease in the body's immune system that mistakenly attacks the joints,[1-5] and cancer.[6-9] ADA inhibitors may change the concentration of adenosine specifically and may act as an anti-inflammatory, anti-HIV and anti-cancer drug, however, usually they have several side effects.[10]ADA inhibitors are classified into two types: transition- and ground-state inhibitors.[11,12] Transition-state inhibitors of ADA resemble a tetrahedral transition-state intermediate as a result of deamination.[11] Ground-state inhibitors of ADA are similar to adenosine. Furthermore, these two types of inhibitors are met with a number of problems, such as toxicity, poor oral absorption, and rapid metabolism.[13,14] ADA has two distinct conformations: closed and open forms, which bind to two types of inhibitors.[15,16] The structural difference is due to an inhibitor-induced conformation change (Fig. 1a and b).
Fig. 1
(a) Active residues in the ADA-PRH complex (PDB Id 1KRM). (b) Active residues in the ADA-FRK complex (PDB Id 1WXY).
Biological macromolecular channels are closely related to the biological function and structural stability of macromolecules and have important structural significance.[17] In some enzymes the active site is relatively exposed to the surface, while in others the active site is deep in the core of the protein. For enzymes that hide the active site, the substrate must pass through a series of tunnels to bind to the active site, a structure that allows more potential protein–ligand interactions to occur than proteins that surface expose the active site. Structurally, an underground active site accessed by one or more tunnels increases the complexity of ligand binding or dissociation. Similarly, the dissociation of inhibitors of enzymes increases their complexity, and through the analysis of protein tunnels, we hope to obtain more effective and safe inhibitors. The active site residues and catalytic mechanism of ADA have been well studied, but the number and characteristics of the tunnels connecting the active site to the protein surface have not been widely characterized.In the early days, drug development, selection, and optimization were mainly based on equilibrium thermodynamic constants such as KD or IC50 values, which to some extent demonstrated the affinity between compounds and the target, but failed to correctly evaluate the in vivo efficacy of the compounds, resulting in the failure of many compound development in phase II clinical trials. Residence time (RT) is an important parameter related to the safety and efficacy of drugs in vivo.[18] The best drug candidate is determined by determining or adjusting the kinetic rate associated with RT.Over the last few decades, computer-aided drug design has emerged to be a powerful technique playing a crucial role in the development of new drug molecules with less time and money.[19] Herein, computer-aided drug design can be used to explore the interactions between ADA and different inhibitors, where it is a useful way to get new excellent inhibitors with low toxicity to regulate the function of ADA.In this study, we used the RAMD method and SMD simulations to explore the inhibitor unbinding tunnels in ADA and each tunnel's crucial residues. Our theoretical results may suggest a clue for ADA investigation.
Materials and methods
Protein preparation
We used two kinds of inhibitors combined with ADA, namely, 6-hydroxy-1,6-dihydro purine nucleoside (PRH) and N-[4,5-bis(4-hydroxyphenyl)-1,3-thiazol-2-yl]hexanamide (FRK). The initial structure of complex in the RAMD simulations was obtained from the RCSB Protein Data Bank (PDB Id 1KRM and 1WXY)[16,20] (Fig. 1a and b). In 1KRM, ADA is presented in a closed form combined with the inhibitor PRH.[20] In 1WXY, ADA is presented in an open form combined with the inhibitor FRK.[16]
Random acceleration molecular dynamics (RAMD) and reaction coordinates
RAMD was created to simplify and facilitate the simulation process and improve the accuracy and reliability of simulation results.[21,22] RAMD is an excellent approach for determining the possible tunnels of a protein, and this technique has been applied to identify new tunnels in many systems.[22-26] RAMD is improved on the basis of conventional simulation and steered molecular dynamics simulation. It applies force in either direction of the ligand to see whether the ligand can be shifted under the action of force. The force exerted on the ligand is expressed as:where k is a constant; r0 is a unit vector to any direction.where m is the number of steps, each force applied to the ligand molecule will maintain m steps, vmin is the minimum ligand velocity under this force, Δt is the time required for a step, rmin represents the minimum displacement of ligand motion in this direction. After m steps, if the average velocity of ligand movement is less than vmin the direction of force will change, and then it will move in a new direction.Two complexes were used in NAnoscale Molecular Dynamics (NAMD) 2.10b1 version software using the CHARMM36 all-hydrogen force field.[27-29] In our initial simulations, m = 10, 20, 30, 40; acceleration = 0.005, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1; and rmin = 0.0005, 0.001, 0.002, 0.005, 0.008, 0.01, 0.02. In order to get better parameters, every combination of the parameters was simulated for 1 ns, and we chose the combination in which the inhibitor travelled in the protein for more than 30 ps. Eventually, we used the following parameters: m = 20; rmin = 0.005; and acceleration = 0.05 successfully ran 300 simulations for each inhibitor. All the results were visualized with Visual Molecular Dynamics (VMD) version 1.9.1.
Steered molecular dynamics (SMD) simulation
SMD is a simulation method for studying ligand (substrate and inhibitor) binding and separation, as well as some mechanical and energy properties in the simulation process.[29] In the simulation of SMD, a time-dependent external force is applied to the ligand, exerting force on the ligand in the set direction, the ligand is finally separated from the receptor. The mechanism of action between ligand and receptor channels can be accurately judged by the time of action, the magnitude of the force, the energy required and the work done. Formula is as follows:where the k is the elastic constants, x⃑ is controlled atomic displacement relative to its initial state, t is the time the action of forces.The two complexes were used in NAMD 2.10b1 version software by using the CHARMM36 all-atom force field.[27-29] Constant-velocity SMD simulations were performed in the present simulations. The pulling velocity was set to 0.01 Å·ps−1. A spring constant of 0.2 kcal mol−1 Å−2 was applied to the inhibitor's centre of mass. SMD simulations were performed from the snapshot structure at 10 ns. To obtain general results and eliminate contingency, we repeatedly stretched each tunnel for three times. The optimal tensile track in parallel experiments in each tunnel was selected for subsequent analysis of tensile force and tunnel residues, and we computed the maximum force for the inhibitor through the unbinding pathway.
Potential of mean force (PMF) with the adaptive-biasing-force (ABF) method
Potential of mean force (PMF) can characterize the substrate's free energy along it, which dictates the speed of transport and the selectivity. Park et al. have proposed a PMF reconstruction method from a limited number of SMD simulation trajectories for a channel system, which used the cumulated expansion to construct the PMF in this study. Formula is as follows:[30,31]where σw2 = 〈W2〉 − 〈W〉2, kB is the Boltzmann constant, T is the temperature. In principle, PMF depth, ΔGoff is inversely proportional to koff and that is to say it is proportional to residence time (RT). The dissociation rate constant (koff) is the constant of the dissociation rate between the reaction ligand and the protein, and the reciprocal residual time (RT) is an important kinetic parameter that directly characterizes the interaction time between the drug molecule and the target protein, which is closely related to the efficacy of the drug. In general, a drug candidate may be more promising if it has a longer RT.Adaptive biasing forces (ABF) can generate quasi-equilibrium trajectories from which the PMF can be deduced. Comparing with other several classic method like umbrella sampling (US), ABF has obvious advantages: ABF introduces the concept of bias force, greatly improves the sampling efficiency, and we don't need to estimate the shape of the free energy barrier in advance and height, also don't need to overlap between adjacent Windows. In this way, the error caused by weight addition and so on will not appear.We choose the atom distance which also used in SMD to complete the ABF-PMF simulations. To enhance the efficiency of the ABF algorithm, the span of the reaction coordinate will be subdivided into equally spaced windows.
Results and discussions
Multiple unbinding pathways identified by RAMD simulations
RAMD, which is an effective approach for identifying the possible tunnel of macromolecules containing protein, has also been successfully applied in many systems.[22-26,32] This study identified the unbinding pathways of two inhibitors (PRH and FRK) from ADA. For the two complexes, 300 simulations were performed to determine the trajectory and observe the unbinding pathway of the inhibitor with VMD. Table 1 shows the number of times each tunnel appears. Three main tunnels (T1–T3) were found in the ADA-inhibitor complex. These three tunnels most likely serve as the unbinding tunnel for inhibitors. Among them, T3 had higher occurrence frequencies for inhibitors unbinding. According to the basic principle of RAMD, we calculated the percentage of occurrences of each tunnel. The unbinding frequencies of T1–T3 in the ADA-PRH system were 17.3%, 17.3%, and 65.3%, respectively, whereas those in the ADA-FRK system were 21.3%, 23.0%, and 55.7%, respectively (Fig. 2a and b). In the ADA-inhibitor systems, T3 had the highest unbinding ratio both in the ADA-PRH and APH-FRK complex. Hence, T3 was the most favourable tunnel for PRH and FRK escape. T2 and T1 were few observed in the simulations.
Summary of the RAMD simulations of ADA-PRH and ADA-FRK
Systems
T1
T2
T3
ADA-PRH
52
52
196
ADA-FRK
64
69
167
Fig. 2
(a) Inhibitor unbinding pathways from ADA as revealed by RAMD simulations. (b) Egress frequency via each tunnel; the upper panel is ADA-PRH, and the lower panel is ADA-FRK.
Tunnel properties on the basis of SMD simulations
According to the result of RAMD, we can see the direction of each tunnel. In order to determine the tunnels direction, the vector direction was determined by taking a point on each protein and center of ligand. After repeated stretching direction confirmation, we finally determined the exact direction of the stretch. In the ADA-PRH complex, the traction direction of T1–T3 is set at the initial location of the ligand in the active site and the centre of the Cα atoms of the residues, Leu137, His17 and Val 16, respectively. In ADA-FRK complex, the directions were set at the inhibitor in the active site and the centre of the Cα atoms of the residues, Tyr102, His17, and Asp296, respectively. We conducted SMD simulations and acquired the trajectory file. We loaded the file into VMD and obtained the force, energy, and distance curve. The rupture forces of the three most frequently observed unbinding tunnels (T1, T2, and T3) were quantitatively estimated. SMD simulations were performed three times for each inhibitor via T1, T2, and T3 (Fig. S1†). SMD was used to further sample both inhibitors and residues along the same unbinding tunnel.[32-37] The maximum force value (Fmax) extracted from the SMD simulation trajectories are listed in Table 2 and Fig. 3a, b. As shown in Fig. 3a and b.
F
max of inhibitors in ADA-PRH and ADA-FRK stretched along three tunnels respectively
Systems
Tunnel
Fmax(pN)
ADA-PRH
T1
9653
T2
10 496
T3
3147
ADA-FRK
T1
2387
T2
2187
T3
2095
Fig. 3
Tension changes of inhibitors leaving each tunnel; (a) ADA-PRH (T1 (black), T2 (red), T3 (blue)). (b) ADA-FRK (T1 (black), T2 (red), T3 (blue)).
The dynamic process of SMD simulation and the change of interaction energy in PRH and FRK dissociation process were also demonstrated and recorded under the best simulation conditions mentioned above. The Lennard-Jones–Short Range (LJ–SR) and the Coulombic Potential–Short Range (Coul–SR) energies from 0 to 10 ns were being analysed to get the electrostatic energy and vdW energy (Fig. 4a–f and S2†). In three tunnels, during the period that the energy starts to change until the time that the all energies go to zero. The higher energy shown in the figure indicated that the bond is weaker and was more likely to disengage. For the ADA-PRH complex, T3 had the highest energy in the three tunnels during 1–2 ns. For the ADA-FRK complex, the energy of T3 was also highest in the three tunnels most of the time during 1.4–8 ns. This result indicated that the inhibitor is more easily released from the T3 tunnel.
Fig. 4
Changes in the interaction energy between ADA and inhibitors during SMD simulation. In ADA-PRH. (a) Energy changes in T1; (b) energy changes in T2; (c) energy changes in T3. In ADA-FRK (d) energy changes in T1; (e) energy changes in T2; (f) energy changes in T3.
PMF calculations
The calculated PMF profiles for ADA-PRH and ADA-FRK complexes were depicted in Fig. 5. The free energy curve revealed the information about unbinding of the two ligands. With the departure of the inhibitor from the initial equilibrium position, the free energy value rapidly increased. In ADA-PRH, the PMF trends of the three tunnels correspond to Fig. 5a–c. The free energy barrier (the PMF depth, ΔGPMF,lowest − ΔGPMF,highest[31]) of the three tunnels were 64.12 kcal mol−1, 65.84 kcal mol−1, 30.54 kcal mol−1, respectively for T1, T2, T3. It's obvious that T3 required lowest energy. The same performance also can be gotten in ADA-FRK complex (Fig. 5d–f), the free energy barrier for T1, T2, T3 were 62.95 kcal mol−1, 66.9 kcal mol−1, 27.88 kcal mol−1. Another way to look at it, the PMF depth ∝ 1/koff and 1/koff was defined as the residence time (RT). Ligands were more likely to escape from T3, and FRK had shorter RT than PRH.
Fig. 5
PMF calculation along each tunnels for PRH (a–c) and FRK (d–f).
Key residues of PRH unbinding pathway along T3, T2, and T1
For the ADA-inhibitor complex, SMD simulations were used to identify the interactions spanning from hydrogen bonds that stabilize the secondary structure or electrostatic interactions between charged side chains at the atomic level to a long-range backbone or domain–domain specific structural properties at the macroscopic level. Once the peak force was obtained, the inhibitors interacted with some residues that may play an important role in inhibitor unbinding pathway. These residues would provide information for designing drugs to easily access the active binding site and enhance or inhibit the activity of the targeted protein.In initial structure of ADA-PRH, PRH was held tightly in the binding pocket by His214, His238, Leu58, His15, His17, and Asp19 with hydrogen bonds (Fig. S3†). Figure shows that His17 made a π–π interaction with PRH.T3 is the favourable pathway for PRH to escape from ADA. At 2.5 ns, the peak point, the hydrogen bonds between Asp19, Asp295, and His238 and PRH were retained, while other hydrogen bonds were broken (Fig. S4†). The π–π interaction between His17 and PRH was broken (Fig. 6a and d). These changes may weaken the interactions between PRH and ADA. However, strong hydrogen bonds with Glu217 and Gly184 were formed with PRH released from ADA. Thus, Gly184, Asp295, His238, and Glu217 were the major residues for the T3 unbinding pathway.
Fig. 6
(a) π–π interactions between His17 and PRH. (b–d) Distance between the centre of the imidazole group of PRH and the centre of the imidazole group of His17 via (b) T1, (c) T2, and (d) T3.
The peak of PRH via T2 unbinding pathway was at 8.5 ns (Fig. 3). All hydrogen bonds disappeared (Fig. S5†). In particular, Glu186 and Glu217 made new hydrogen bonds with PRH. The π–π interaction between His17 and PRH was broken (Fig. 6a and c). Gly216 formed one hydrogen bond with PRH. These hydrogen bonds interactions may prevent the release of PRH escaping from ADA. Thus, the major residues in T2 were Glu186, Glu217, and Gly216.The active residues for PRH at the peak time (∼7.5 ns) in T1 were shown in Fig. S6.† At this time, all hydrogen bonds at the active pocket disappeared. The π–π interaction between His17 and PRH was broken (Fig. 6a and b). Moreover, new hydrogen bonds were formed (Glu186 and Arg156 made two hydrogen bonds with PRH, and Glu217, Val218, Leu58, and Asp185 formed a hydrogen bond with PRH). These hydrogen bonds may prevent the release of PRH from ADA via T1. Hence, Glu186, Arg156, Glu217, Val218, Leu58, and Asp185 were the important residues in T1.In particular, His17 made a π–π interaction with PRH in the active pocket, and this interaction disappeared when PRH escaped from ADA via T3, T2, and T1 (Fig. 6a–d). Hence, PRH escaped from ADA via T3 and may break through less obstructions than those via T2 and T1.
Key residues of FRK unbinding along T3, T2, and T1
FRK was held in the binding pocket by Asp19 with a hydrogen bond initially (Fig. S7†). At the peak time (∼1.9 ns), the hydrogen bond between Asp19 and FRK was retained (Fig. S8†). In particular, Trp117 made a strong π–π interaction with FRK (Fig. 7a). These interactions (hydrogen bond and π–π interaction) may prevent the release of FRK leaving from ADA. Thus, Asp19 and Trp117 were the major residues for FRK via the T3 unbinding pathway.
Fig. 7
(a) π–π interactions between Trp117 and FRK in T3. (b) σ–π interactions between Leu62 and FRK in T2. (c) π–π interactions between Phe65 and FRK in T2. (d) σ–π interactions between Leu62 and FRK in T1.
As shown in Fig. 3b, the peak of FRK via the T2 unbinding pathway was at 1.9 ns. The hydrogen bond between Asp19 and FRK was retained (Fig. S9†). In particular, Phe65 formed a π–π interaction with FRK (Fig. 7c). Leu62 made an σ–π interaction with FRK (Fig. 7b). These interactions may prevent the release of PRH escaping from ADA. In T2, the major residues were Asp19, Phe65, and Leu62. The active residues for FRK via T1 were observed (Figs. S10† and 3b) at the peak time (∼2.0 ns). At this time, the hydrogen bond between Asp19 and FRK at the active pocket disappeared. Two new hydrogen bonds between Glu217 and Asp296 were formed. Leu62 made an σ–π interaction with FRK (Fig. 7d). These interactions may prevent the release of FRK from ADA via T1. Thus, Glu217, Asp296, and Leu62 were the important residues of FRK escaping from ADA via T1.A strobe image of PRH escaping from ADA via T3, T2, and T1 tunnel was shown Fig. 8a–c. According to the above results, PRH may have escaped from ADA via T3 and may break through less obstructions than those via T2 and T1.
Fig. 8
Strobe images of PRH (a–c) and FRK (d–f) escaping from ADA via T3, T2, and T1.
FRK escaped from ADA via T3, T2, and T1 tunnel, as shown in Fig. 8d–f. FRK may have escaped from ADA via T3 and may break through less obstructions than those via T2 and T1.
Alanine scanning for predicting the key residues in the peak of the three tunnels
Alanine scanning mutagenesis was used to analyze the function of specific amino acid residues on the surface of protein. SMD trajectories were employed to carry out binding free energy calculations and perform an alanine scanning of the peak site residues in six complexes by using FoldX software[38] (Fig. 9a–c). The most active residues that remarkably contributed to the binding affinity was Trp117 for PRH-ADA via T3, Phe65 for PRH-ADA via T2, and Leu62 for PRH-ADA via T1. However, for the FRK-ADA at peak state, the results of alanine scanning revealed the energetic hot spots corresponding to different residues. In addition, in terms of same tunnel, different inhibitor in protein, Leu58 in T1, Glu186 in T2, Asp295, Glu257 and Gly184 in T3 showed wide difference. This interesting phenomenon can be further explored.
Fig. 9
Computational alanine scanning of the peak site residues of six complexes. Binding free energies and alanine scanning of residues are shown for (a) T1, (b) T2, (c) T3 in two systems. PRH-ADA is shown in purple bars and FRK-ADA is shown in blue bars.
Conservation of tunnel residues
It was much known that residues conservation rate in the superfamily of the protein is important for studying the interactions between protein and inhibitors. ADA belongs to adenosine/adenosine monophosphate (AMP) deaminase family.[39] This protein family concludes ADA, ADA2, adenosine deaminase-like protein (ADAL), AMP deaminase (AMPD1, AMPD2 and AMPD3).[40] The complete protein sequences of each protein were obtained from UniProt Knowledgebase and aligned via BioEdit software (version 7.2.5).[41] The complete amino acid sequence alignment diagram is shown in Fig. S11† and the shade threshold was set as 50%. In the above, we got the residues in the three tunnels that had significant effect on inhibitor binding. To explore their conservatism in the family, we extracted fragments of each residue for comparison and calculated the identity percentage,[42] as shown in the Table 3. The identity percentage of each amino acid in ADA is calculated by adding 1 to the number of proteins with the same residue as ADA and dividing by the total number of proteins. It can be concluded that Trp117, Glu186, His238, Asp295 are conservative in adenosine and AMP deaminases family proteins.[8,43,44] Hence, Trp117 for PRH-ADA via T3 may play important roles in the unbinding tunnel for inhibitor leaving.
Retention of ADA channel residues in other proteins
ADA
ADA2
ADAL
AMPD1
AMPD2
AMPD3
Percent identity
Asp19
Cys157
Leu75
Ala316
Ser369
—
0
Leu58
Tyr219
Leu138
Phe378
Phe431
Phe60
33%
Leu62
Asp223
Leu142
Asp382
Asp435
Asp64
33%
Phe65
Trp226
Phe145
Asn385
Asn438
Asn67
33%
Trp117
Trp286
Trp197
Trp445
Trp498
Trp127
100%
Arg156
Leu322
Arg236
Val491
Ile544
Ile173
33%
Gly184
Gly348
Gly264
Asp522
Asp575
Asp204
50%
Asp185
Arg349
Asp265
Asp253
Asp576
Asp205
83%
Glu186
Glu350
Glu266
Glu524
Glu577
Glu206
100%
Gly216
—
Gly296
Pro580
Pro633
Pro262
33%
Glu217
His378
Glu297
His581
His634
His263
33%
His238
His406
His318
His603
His656
His285
100%
Asp295
Asp463
Asp375
Asp658
Asp711
Asp340
100%
Conclusions
ADA is involved in numerous purine metabolism processes. In previous studies, the unbinding tunnel of two different inhibitors remains unknown. In this study, we used RAMD and SMD simulations to explore the unbinding tunnel of two inhibitors (PRH and FRK) via ADA. Three main tunnels were most frequently observed, namely, T1, T2, and T3. Among them, T3 was the most favourable unbinding tunnel for two inhibitors. We also predicted and compared the critical residues in T1, T2, and T3 to ADA. Our results provide a reference for further studies on ADA and a theoretical basis for the development of subsequent ADA inhibitors.
Authors: Doris A Schuetz; Wilhelmus Egbertus Arnout de Witte; Yin Cheong Wong; Bernhard Knasmueller; Lars Richter; Daria B Kokh; S Kashif Sadiq; Reggie Bosma; Indira Nederpelt; Laura H Heitman; Elena Segala; Marta Amaral; Dong Guo; Dorothee Andres; Victoria Georgi; Leigh A Stoddart; Steve Hill; Robert M Cooke; Chris De Graaf; Rob Leurs; Matthias Frech; Rebecca C Wade; Elizabeth Cunera Maria de Lange; Adriaan P IJzerman; Anke Müller-Fahrnow; Gerhard F Ecker Journal: Drug Discov Today Date: 2017-04-13 Impact factor: 7.851
Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9