Long-Timescale Simulations Revealed Critical Non-Conserved Residues of Phosphodiesterases Affecting Selectivity of BAY60-7550.

A major obstacle of the selective inhibitor design for specific human phosphodiesterase (PDE) is that highly conserved catalytic pockets are difficult to be distinguished by inhibitor molecules. To overcome this, a feasible path is to understand the molecular determinants underlying the selectivity of current inhibitors. BAY60-7550 (BAY for short; IC50 = 4.7 nM) is a highly selective inhibitor targeting PDE2A which is a dual-specificity PDE and an attractive target for therapeutic intervention of the central nervous system (CNS) disorders. Recent studies suggest that molecular determinants may be in binding processes of BAY. However, a detailed understanding of these processes are still lacking. To explore these processes, High-Throughput Molecular Dynamics (HTMD) simulations were performed to reproduce the spontaneous association of BAY with catalytic pockets of 4 PDE isoforms; Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD) simulations were performed to reproduce the unbinding-rebinding processes of FKG and MC2, two pyrazolopyrimidinone PDE2A selective inhibitors, in the PDE2A system. The produced molecular trajectories were analyzed by the Markov state model (MSM) and the molecular mechanics/generalized Born surface area (MM/GBSA). The results showed that the non-covalent interactions between the non-conserved residues and BAY, especially the hydrogen bonds, determined the unique binding pathways of BAY on the surface of PDE2A. These pathways were different from those of BAY on the surface of the other three PDE isoforms and the binding pathways of the other two PDE2A inhibitors in PDE2A systems. These differences were ultimately reflected in the high selectivity of this inhibitor for PDE2A. As a result, this study demonstrates the critical role of the binding processes in the selectivity of BAY, and also identifies the key non-conserved residues affecting the binding processes of BAY. Thus, this study provides a new perspective and data support for the further development of BAY-derived inhibitors targeting PDE2A.

Introduction

Human phosphodiesterases (PDEs) comprise of 11 enzyme families, namely PDE1 ∼ PDE11. They can degrade the phosphodiester bonds of cAMP or cGMP to control the concentrations of these second messengers [1]. As second messengers are critical intracellular signaling molecules, the dysfunction of PDEs can lead to a variety of serious diseases, such as cancer [2], heart failure [3], and schizophrenia [4]. As a result, PDEs are potential therapeutic targets and have been investigated for years.

The high selectivity of inhibitors means low off-target rates, resulting in fewer side effects. Since PDE isoforms involve different physiological activities in various human tissues, the selectivity is an essential property of PDE inhibitors. The homology of PDE catalytic domains is ranging from 20 % to 45 %, suggesting the possibility to design selective inhibitors target specific isoforms [5]. At present, there are selective inhibitors that can be used in clinical applications targeting PDE3, PDE4, PDE5, and PDE10 [6], [7], [8], [9], [10], [11], [12]. However, they all have more or less side effects. These include the increased mortality, headache, and vision loss caused by Milrinone (PDE3 inhibitor), Roflumilast (PDE4 inhibitor) and Sildenafil (PDE5 inhibitor), respectively [6], [7], [8], [9], [10], [11]. One possible reason is that the catalytic pockets of PDE isoforms are highly similar in structure and sequence, which makes these pockets difficult to be distinguished by inhibitors (Fig. 1A and S1). So, there are still urgent medical needs to develop highly selective inhibitors targeting specific PDE isoforms. In order to improve the selectivity of inhibitors, a feasible path is to find the molecular determinants of the selectivity of currently available inhibitors and make full use of them.

Fig. 1

(A) Conservation of residues of 4 PDEs (PDE2A, PDE4B, PDE5, PDE10A) according to the Clustal sequences alignment (See Fig. S1). (B) Structural formulae of BAY, FKG and MC2. (C) Initial structure of PDE-BAY complex for HTMD.

BAY60-7550 (BAY for short; IC50 = 4.7 nM) is a highly-selective inhibitor targeting PDE2A (Fig. 1B) [13]. Its selectivigy targeting PDE2 is 50 times that of PDE1C and more than 100 times that of other PDE subfamilies. Over the past decade, researchers have taken advantage of the selective inhibition of PDE2A by BAY to delve into the function of PDE2A. Studies have shown that BAY can significantly increase the concentrations of cAMP and cGMP in stimulated primary neuronal cultures, revealing the dual hydrolytic function of PDE2A [13]. In the mouse behavior model, BAY could significantly improve the cognition and communication ability of rats, indicating that PDE2A is related to the memory and learning ability of mammals [14]. BAY showed an antidepressant effect in anxious mice, indicating that PDE2A can be used as a drug target for psychiatric disorders [15].

Although BAY has not entered clinical use due to its moderate pharmacokinetic properties [14], the molecular mechanism of its high selectivity is worthy of investigation for the development of novel PDE2A inhibitors. For a long time, the affinity has been an indicator to characterize the inhibitor selectivity [16]. In molecular docking studies, this molecule showed similar conformation and affinity in the catalytic pockets of all PDE isoforms [17]. This implied that the selectivity of BAY may not entirely dependent on its affinity. In pharmacology, affinity is usually defined as the concentration of ligand required to occupy 50 % of the targets at equilibrium [18]. In fact, the in vivo concentration of inhibitor molecules changes at a rate much greater than their binding rate. As a result, the above equilibrium is difficult to achieve in the body [16]. Therefore, it is not enough to evaluate the selectivity of BAY only by its affinity.

Moreover, some studies have shown that the kinetic properties of inhibitors, such as the residence time, also affected their selectivity [16], [19]. On the other hand, in our previous in silico study, this inhibitor took various binding pathways on the surface of PDE isoforms, resulting in different binding poses in the pockets of 6 PDEs; moreover, along the binding pathway in the PDE2A system, BAY can enter the catalytic site with the highest probability and form the simulated pose approximating the crystal one [20]. Therefore, the binding pathways of BAY on the surface of PDEs should be an important factor in determining the selectivity of this inhibitor.

In the absence of external forces, the binding pathways of BAY should be mainly determined by the interactions between the inhibitor and the surface residues outside catalytic sites of PDEs. For any PDE isoform, the active pocket mainly consists of conserved and similar residues, while the region outside the pocket contains a large number of non-conserved residues (Fig. 1A) [21]. Compared with conserved and similar residues, non-conserved residues have significantly different side chains. So, the non-covalent interactions between BAY and non-conserved residues should be the reason for the different behaviors of the inhibitor on the surface of PDEs. These behavioral differences may lead BAY to choose different binding pathways. Therefore, the non-conserved residues, especially those outside catalytic sites, may be critical for the binding pathways of BAY.

However, there is still a need to achieve a detailed knowledge of the interactions between BAY and non-conserved residues of PDEs. Molecular dynamics (MD) simulation is an effective method to study the binding processes of drug molecules. This method describes the time-dependent motions of the drug-protein systems at the atomic level, which can reproduce the binding or unbinding processes of drugs [22], [23]. According to IC50 values, the selectivity of BAY for PDE2A is about 389, 150 and 200 times higher than those for PDE4B, PDE5 and PDE10A, respectively [13]. As mentioned above, these four PDEs have been identified as targets for the treatment of the corresponding diseases, and some selective inhibitors targeting them have also been in the clinical use. Moreover, the crystal structures of their catalytic domains are relatively complete. Therefore, these four PDE isoforms were selected for comparative studies of the interactions between non-conserved residues and BAY. In this study, High-Throughput Molecular Dynamics (HTMD) simulations at microsecond scale were used to reproduce the association processes of BAY on the surface of those four PDE isoforms [24]. In HTMD, BAY started at a point far away from the pockets of PDEs and spontaneously bound to these pockets driven by its interactions with the residues of PDEs (Fig. 1C). In order to avoid the well-known sampling problem of long-term trajectories, 10 echoes of 25 ns unbiased molecular dynamics simulations were carried out in an iterative stepwise manner for each PDE-BAY system (Fig. 2) [25], [26]. Only in the PDE2A-BAY system, BAY successfully bound into the catalytic site and formed the simulated pose approximating the crystal one (RMSD = 3.51 Å). All binding processes were characterized by means of the Markov state model (MSM) [27], [28] and the molecular mechanics/generalized Born surface area (MM/GBSA) [29], [30]. The data of MSM and MM/GBSA showed that (1) the non-conserved residues outside the pocket may affect the conformational preference of BAY in the binding processes through the hydrogen bonds with BAY; (2) hydrogen bonds between E846 and BAY nay directly determin that BAY adopts a crystal-like pose to bind in the catalytic site PDE2A to a 78.1 % degree. (3) the metoxyphenyl and pyrazolopyrimidinone groups of BAY were involved in 3.5 % and 74.6 % of the hydrogen bonds with E846, respectively; (4) by interacting with the non-conserved residues outside the pocket of PDE2A, the metoxyphenyl group can promote hydrogen bonds between E846 and the pyrazolopyrimidinone group. Thus, the hydrogen bonds between E846 and BAY may facilitate the successful binding of this inhibitor to the PDE2A catalytic site and induce BAY to adopt a simulated binding pose that approximates the crystal one.

Fig. 2

Workflow of HTMD, LiGaMD and data analyses.

Given the critical hydrogen bonds between the pyrazolopyrimidinone group and E846, an interesting hypothesis is that this hydrogen bond may also be present in the binding processes of other selective pyrazolopyrimidinone PDE2A inhibitors. Compared with the unbiased MD, LiGaMD has more efficient sampling efficiency of inhibitor binding pathways by applying specific potential boosts to the complex and inhibitor, respectively. So, to test the above hypothesis, the microsecond-scale Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD) simulations were used to reproduce the unbinding and rebinding processes of FKG (IC50 = 2 nM) and MC2 (IC50 = 1.3 nM), the other two pyrazolopyrimidinone inhibitors targeting PDE2A (Fig. 1B) [5], [31]. Probably due to the lack of the metoxyphenyl group, FKG and MC2 hardly formed hydrogen bonds with non-conserved residues during their binding processes, thus leading to their different binding pathways from those of BAY.

In summary, this study illustrates the critical role of the binding processes for BAY selectivity, and also identifies non-covalent interactions between non-conserved residues and inhibitors that can determine the BAY binding poses. In this way, this study provides a new perspective and the data support for the further development of BAY-derived inhibitors targeting PDE2A.

Materials and methods

Construction of PDE-BAY complex systems for HTMD

The catalytic domains of PDE2A, PDE4B, PDE5, and PDE10A (PDB ID: 4HTX, 4KP6, 2H42, 2OUP) were used to construct the HTMD simulation systems [17], [32], [33], [34]. First of all, the centroid of each PDE isoforms was located at the origin of coordinates, and the entrance of the catalytic pocket faced the positive direction of Y-axis. To ensure adequate conformational adjustment before entering the catalytic pockets, BAY60-7550 was put at least 45 Å and 9 Å from the catalytic site and the PDE surface, respectively (Fig. 1C). The PDE-inhibitor complex was immersed in a rectangular water box containing 150 mM NaCl. PDE, BAY60-7550 and water molecules were described by the CHARMM36 Force Field, CHARMM General Force Field (CGenFF) and TIP3P water model, respectively [35], [36], [37].

Our previous study has shown that BAY60-7550 enters the catalytic pocket of PDE2A in the first 25 ns at 320 K, and the binding pathway hardly passes the PDE surface located at the back of the pocket [20]. Therefore, 25 ns was set as the time of the short-time scale MD simulation, and the residues on the back of the pocket were constrained to the initial coordinates by the harmonic constraints. This reduced unnecessary translations and rotations of PDEs. On the other hand, the motion range of BAY60-7550 was confined to a cubic box centered on the pocket entrance with an edge length of 50 Å (Fig. 1C). This increased the probability of BAY60-7550 finding the pocket entrance. Under periodic boundary conditions, the MD simulation of the above five systems were carried out by ACEMD (Ver.2.0) [38]. The temperature was controlled by the Langevin Thermostat at 320 K in the NVT ensemble [39]. The simulation step size was 4 fs. The simulation frames were output very 100 ps. The Particle Mesh Ewald (PME) method was used to deal with electrostatic forces at long distances [40].

Workflow of HTMD

In this study, 10 echoes of MD simulations were performed using HTMD. Each echo produced 10 to 15 short-time scale trajectories. These trajectories were used to build the MSM of this echo and generate the start poses of the next echo according to the adaptive sampling algorithm. Previous studies have shown that: 1) the metal ions in its catalytic pocket are critical for the catalytic activity of PDEs [41], [42]; 2) the hydrophobic subpocket (H-pocket) and the π-π stacking with F862 (according to the serial number of PDE2A) are important for the inhibitor affinity [17].Therefore, the metrics for building the MSM were set to be the distances between atoms of PDEs (includes these metal ions and F862 side chain heavy atoms) and four specific heavy atoms of BAY60-7550. The above operations were automatically completed by HTMD (Ver.1.0) in an iterative stepwise manner. For each PDE-inhibitor system, 10 echoes of MD simulations yielded at least 100 short-time scale trajectories, with a total simulation time of at least 2.5 μs (Fig. 2).

Construction of PDE2A-FKG/MC2 complex systems for LiGaMD

FKG and MC2 were named in this study according to the contents of the corresponding PDB files. The catalytic domains of PDE2A (PDB ID: 6CYD, 4HTX) in complex with FKG or BAY were used to construct the PDE2A-inhibitor system for the LiGaMD simulation [5]. Since the complex structure of MC2 and PDE2A was not included in the Protein Data Bank, the docking pose of MC2, which was the most crystal-like conformation, was used to construct the PDE2A-MC2 system together with the catalytic domain structure of PDE2A (PDB ID: 4HTX). Each of these complexes was immersed in a rectangular water box containing 150 mM NaCl. PDE2A, inhibitor and water molecules were described by the AMBER14SB Force Field, Generalized Amber Force Field (GAFF), and TIP3P water model, respectively [36], [43], [44]. These systems were all energetically minimized for 2000 steps and equilibrated at 310 K and 1 atm pressure for 12 ns. Then, 2 ns unbiased MD was applied to each of these systems to calculate LiGaMD acceleration parameters. After 50 ns LiGaMD equilibration, at least 6 independent 500 ∼ 2000 ns LiGaMD production simulations with randomized initial atomic velocities were performed on the PDE2A-inhibitor systems. All LiGaMD simulations were run at the “dual-boost” level. When the distance between heavy atoms of inhibitors and the side chain of F862 was more than 4 Å (distance in crystal pose), the boost potential was applied to the ligand non-covalent interaction potential energy. The remaining potential energy of the entire system was boosted constantly. The simulation frames were saved every 100 ps for analysis. For both systems, LiGaMD simulations yielded trajectories with a total simulation time of at least 7.7 μs. LiGaMD simulations were carried out by the AMBER18 software package [45].

Discretization of MD trajectories

As described below, the trajectories generated by HTMD and LiGaMD would be used to construct MSMs to obtain binding pathways and representative poses of inhibitors. MSM is a dicrete-state stochastic model, while the MD trajectories is continuous. Therefore, reasonable discretization of MD data is needed to balance the statistical error and the system approximation error when constructing MSM [26]. Despite many efforts, there is still no general method to discretize high-dimensional data. Thus, according to data characteristics and scientific problems, the following methods were adopted in this study.

First, after removing the solvent and periodic boundary conditions, the PDE-inhibitor conformation in all frames was superimposed onto the crystal conformation of the PDE-inhibitor complex. Taking the crystal pose as the reference, the smaller RMSD of the simulated pose of inhibitor was, the closer the simulated pose was to the crystal one. Since this study focused on the interactions between PDE and inhibitors in bound states, frames containing free-state inhibitors were removed from the MD trajectories. The remaining trajectories were then used for the k-means clustering based on the heavy atom coordinates of inhibitors. In the clustering process, the number of categories kept increasing until the average radius of categories was less than 2 Å. Lone frames and frames belonging to categories with radius greater than 2 Å were removed from MD trajectories. In this way, the processed trajectories not only retain the original time order of all MD frames, but also ensured that the poses of inhibitors belonging to the same categories had a high similarity. This was beneficial to reduce the statistical errors of metastates and the system approximation errors of the inhibitor binding processes in the subsequent MSM construction. Discretization of MD trajectories were carried out by R-3.3.4 and VMD software package [46], [47].

Construction and analysis of MSM

The discrete trajectories of PDE-inhibitor complexes were used to construct the MSMs carried out by the PyEMMA-2.5.7 software package [48]. The metrics of MSMs were the heavy atom spatial coordinates of inhibitors. After the dimension reduction and clustering, the discrete HTMD or LiGaMD trajectories were divided into several metastates. As the coarse representations of MSMs, each metastate contained 10 representative conformations of PDE-inhibitor complexes. In these conformations, the mean RMSD of the inhibitor with respect to its crystal pose was defined as the RMSD of the corresponding metastate. In this way, the smaller the RMSD of a metastate, the closer the inhibitor in that metastate was to the crystal pose. Therefore, the state with the smallest RMSD was the end state. According to the transition path theory (TPT) [49], [50], the pathways from each state to the end state were found, and the fluxes of these pathways were calculated. The flux of a pathway could be understood as the probability that the inhibitor reached the end state along that pathway. On this basis, the flux of a metastate in this study was defined as the probability that this metastate was passed by other metastates among all the pathways towards the end state. For metastate MJ, its flux was calculated by the following formula:Where, M was the total number of metastates, Km was the number of pathways through metastate MJ in these pathways from metastate m (excluding the end state and metastate MJ) to the end state, and Pm was the flux of the Kth pathway in these pathways.

According to the above definition, the initial state was determined to be the metastate with a high RMSD and the minimum flux. The pathways from the initial state to the end state should be all possible pathways for binding or unbinding of inhibitors. Then, in this study, the effective flux of a metastate was defined as the probability that the inhibitor passes through this metastate along the above pathways. The effective flux of state MJ was calculated by the following formula:Where, Ki wass the number of pathways through state MJ in the pathways from the initial state (excluding the end state and state MJ) to the end state, and Pi was the flux of the Kth pathway in these pathways. In this study, for each PDE-inhibitor system, the dominant pathways were defined as those with high fluxes and the sum of their fluxes greater than 60 %. The metastates included in the dominant pathways were used for the binding free energy calculation and binding kinetics analysis of inhibitors.

Binding free energy calculations using MM/GBSA

In order to obtain the contributions of individual amino acids to.

the inhibitor binding energy, MM/GBSA coupled with per residue energy decomposition was used to calculate the binding free energy between PDEs and inhibitors in all metastates [29], [30]. For each metastate, one representative conformation was used to construct the aqueous solution systems of PDE-inhibitor complex. After the equilibrium under the NPT ensembles, the last 10 frames were used for the MM/GBSA calculation. In this way, the ten residues that contribute most to binding energies were identified for each metastate, The conservations, free energy contributions and side chain poses of the top 10 residues were used to identify residues that determine inhibitor selectivity. MM/GBSA was carried out using the AMBER18 software package [45].

Results

Binding pathways of BAY in PDE2A system

For each PDE-BAY system, the HTMD produced at least 100 short-time scale trajectories, with a total simulation time of at least 2.5 μs. In these trajectories, the RMSDs of BAY with respect to its crystal pose dropped below 15 Å for many times (Fig. S2). The depth from entrance to the catalytic site of PDE pockets was known to be 15 Å [51]. With the above criterion, the RMSDs showed that BAY can spontaneously bind into the pocket in four PDE-BAY systems.

Our previous study showed that the RMSD of BAY relative to its crystal pose remained decreasing or stationary for at least 100 ns once this RMSD falled below 5 Å in the unbiased MD simulations. This indicated that BAY was approaching its crystal pose intermittently during this 100 ns. This trend may be difficult to be shown in the 25 ns unbiased MD trajectories of HTMD in this study. In fact, the π-π stackings between the pyrazolopyrimidinone group of three inhibitors and the phenyl ring of F862 side chain are the main features of their crystal pose; and the π-π stacking was also considered to be one of the key factors driving BAY binding into the catalytic site [52]. In other words, this π-π stacking is necessary for BAY binding to the catalytic site. Ji-hua Deng et al. showed that the maximum distance between the two rings of π-π stacking is 3.965 Å [53]. As shown in Fig. S2 and S3, when the RMSDs of inhibitors were below 5 Å, the distances between pyrazolopyrimidinone group and the phenyl ring of F862 side chain often also fell below 3.965 Å. Thus, RMSD less than 5 Å was considered as the criterion to judge whether three inhibitor successfully bound to the catalytic site.

Among the above four systems, only in the PDE2A system did the RMSDs of BAY drop below 5 Å, with a minimum RMSD of 3.51 Å (state 10 in Fig. 4 and Fig. S2). The corresponding simulated pose exhibited the same orientation and essentially the same spatial position as the crystal one. Therefore, only in the PDE2A system, BAY can spontaneously bind into the catalytic site and form the simulated crystal-like poses.

Fig. 4

Nine metastates comprising the top 5 binding pathways that dominate the association processes of BAY in PDE2A system. For each metastate, upper is the Top10 residues with their contribution to BAY binding energy, and bottom is the poses of BAY (violet stick model) and Top10 residues (stick model) of PDE2A (cartoon model). The conservation of residues is indicated by their colors. Hydrogen bonds are displayed as yellow dotted lines. Metastates are displayed in descending order of RMSD from left to right and from top to bottom. Simulated pose (violet stick model) of BAY is aligned to its crystal pose (cyan stick model) in the end state.

As mentioned in the introduction, the spontaneous binding processes of BAY may be determined by the interactions between the inhibitor and non-conserved residues on the surface of PDEs. To verify this hypothesis, the PDE2A system was analyzed as follows. In MSM, the discretized trajectories of PDE2A-BAY complex were divided into 24 metastates, and the state 10 with the minimum RMSD was determined as the end state. The MSM also gave all possible pathways from the other 23 states to state 10 and the corresponding fluxes. On this basis, the effective flux of each metastate was given by Formula (2) in Section 2.5. Among the four metastates with large RMSDs, state 3 with the smallest flux was set as the initial state. From state 3 to state 10, there were 44 possible binding pathways, among which the first 5 pathways whose flux sum was 64.9 % constitute the dominant pathways of the system. That is, BAY had a 64.9 % chance of successfully binding into the catalytic site along these five pathways. Moreover, the metastates making up these dominant pathways participated in all binding pathways mentioned above (Fig. 3C). That is, BAY had 100 % probability of successfully binding to the catalytic site through these 9 metastates. Therefore, the metastates that comprised the dominant pathways should contain information about the molecular determinants of the binding selectivity of BAY.

Fig. 3

(A) Effective flux (left bar), energy contribution from non-conserved residues (right bar) and RMSD (broken line) of 9 metastates comprising the top 5 binding pathways that dominate the association processes of BAY in PDE2A system. (B) Center of mass (COM) of BAY in those 9 metastates. COMs are displayed as spheres and their colors follow the bars in (A). (C) Flux network from the initial state (state 3) to the end state (state 10). Colors of metastates and their fluxes follow the bars in (A). Values of fluxes from state 3 to state 6 and 14 are labeled on the left. Values of fluxes from state 18, 21, 22 and 23 to state 10 are labeled on the right. Non-dominant states are displayed in grey.

There were nine metastates (state3, 6, 14, 22, 21, 23, 18, 24, 10, in descending RMSD order) that made up these dominant pathways (Fig. 3A). The contribution ratios of non-conserved residues to BAY binding free energy were less than 40 % in state 14 and 10, and less than 20 % in the other eight metastates. This implied that the non-conserved residues did not dominate the affinity for BAY. That is, the affinity should not be the main way in which non-conserved residues influence the movement behavior of the inhibitor on the PDE2A surface.

In addition to the affinity, the conformational preference of a molecule should be another critical determinant of its biological properties, which in turn affect its behavior [54]. In order to conduct conformational preference analysis for BAY, as shown in Fig. 4, the top 10 residues contributing to BAY binding energy in each metastate (hereinafter referred to as Top10 residues) were displayed as the stick model around BAY, and their binding energy contribution values were shown in the upper bar chart. In the initial state, BAY was between the alpha helices H10 and H11 (Fig. 3B). This indicated that the region between these two helices was the most likely starting point of the binding pathways of BAY. In this state, non-conserved residues contributed 10.9 % of the binding energy to BAY, but accounted for half of the Top10 residues (Q794, R762, I776, R775, I764; Fig. 4). Among them, Q794, with the second highest free energy contribution, and D763, a weakly similar residue, jointly formed hydrogen bonds with the pyrazolopyrimidinone group of BAY (Fig. 4). Because the hydrogen bond is a short-range, directional and favorable interaction, both ligand and receptor need to make appropriate conformation to achieve the specific bonding by the two hydrogen bonds [55]. In addition, the local electrostatic potential influenced by R762 and R775 with positive side chains may also affect BAY recognizing the binding site in state 3, although they may be less directional than hydrogen bonds. The electrostatic interactions above may require BAY to take the specific pose as shown in the Fig. 4 to bind between H11 and H10 in state 3. In this pose, the phenyl group of BAY was towards the PDE2A pocket entrance. This may facilitate the interaction of this hydrophobic group with hydrophobic residues near the entrance in subsequent states. Thus, in state 3, non-conserved residues not only contributed a small amount of affinity, but more importantly assisted BAY conformational adjustment through electrostatic interactions such as hydrogen bonds.

As the inhibitors in state 6 and 3 were spatially close to each other and their RMSDs were both greater than 20 Å (Fig. 3A), these two metastates were classified as the early states on the binding pathways of BAY in this study. In all pathways state 6 followed state 3. This implied that the two metastates were strongly related in time and space. Non-conserved residues contributed 17.8 % of the binding energy to BAY, but accounted for 7 of the Top10 residues (Fig. 4). Among them, R762, I776, R775, and I764 both appeared in state 3 and 6. Moreover, R762 and the weakly similar residue D763 formed hydrogen bonds with the pyrazolopyrimidinone group of BAY. Thus, similar to state 3, non-conserved residues not only contributed a small amount of affinity, but more importantly assisted BAY conformational adjustment through hydrogen bonds..

Of the flux from the initial state, 53.3 % flowed to state 6, and the rest flowed to state 14 (Fig. 3C). The binding site of BAY in state 14 was between the alpha helices H15 and H12 above the pocket entrance, and the molecule did not enter the pocket (RMSD greater than 15 Å; Fig. 3A). Among the 9 metastates that made up the dominant pathways, except the two early states, all the others were downstream states of state 14 (Fig. 3C). Therefore, this state was classified as the middle state in this study. State 14 had the highest effective flux (65.3 %). This indicated that BAY had a 65.3 % probability of successfully binding to the catalytic site via state 14. Excluding the end state, state 14 had the highest (30.6 %) binding energy contribution ratio of non-conserved residues (Fig. 3A). Seven of the Top10 residues in state 14 were non-conserved. Among them, Q781 formed two hydrogen bonds with the metoxyphenyl group of BAY to assist the conformational adjustment of this inhibitor. Other non-conserved residues were hydrophobic except K778. The hydrophobic interactions of F777 and P869 with metoxyphenyl and phenyl groups contributed to the first two binding free energies per residue (Fig. 4). Guided by these interactions, BAY adopted a pose with the pyrazolopyrimidinone group facing outward, the phenyl group near the entrance of the pocket, and the metoxyphenyl group toward the hydrophobic region around F777. As mentioned above, the BAY poses in the downstream metastates had a 65.3 % probability of being directly or indirectly transformed from this pose. Therefore, in state 14, non-conserved residues contributed 30.6 % of the affinity of PDE2A for BAY, and also assisted BAY to adopt the appropriate conformation through hydrogen bonds and hydrophobic interactions, thus directly and indirectly determining the most of the behaviors of BAY in its downstream metastates.

In state 21, 22, 23, 24, and 18, inhibitors were located at the M−loop, and their RMSDs were between 11 and 15 Å (Fig. 3A and 3B). That is, the BAY poses in these metastates were closer to the crystal pose than those in the early and middle states. So, these metastates were classified as the late states in this study. The total flux of the binding pathways with these late states was up to 83.4 %. This indicated that the late states dominated the binding processes of BAY. According to MSM theory, the pose of BAY in the end state depended only on the next-to-last state, not on the earlier ones. In fact, the 5 late states as the next-to-last state were in most binding pathways, and the total flux of these pathways is also 83.4 %. In other words, there was an 83.4 % probability that the crystal-like pose of the end state is directly converted from the BAY pose in the late states.

As those in the early states, the non-conserved residues in the late states did not dominate the affinity for BAY. However, BAY formed a total of 10 hydrogen bonds with PDE2A in the five late states; half of these hydrogen bonds in the three late states (state 21, 22 and18) were formed between pyrazolopyrimidinone group of BAY and the non-conserved residue E846 (Fig. 4). Moreover, the total flux of the binding pathways with these three late states was as high as 74.8 % (Fig. 3C). In other words, in the process of BAY reaching the end state, the inhibitor had at least a 74.8 % probability of forming hydrogen bonds with E846 at the M−loop by virtue of its pyrazolopyrimidinone group. This meant that the hydrogen bonds between E846 and pyrazolopyrimidinone group were critical to the processes by which BAY bound into the catalytic site of PDE2A and formed the crystal-like pose.

As mentioned above, both ligand and receptor need to adopt the specific conformations in order to form hydrogen bonds for the specific binding. In the late state 22, 21 and 18, the hydrogen bonds between E846 and the pyrazolopyrimidinone group made BAY adopt the poses in which their pyrazolopyrimidinone groups appeared above the M−loop and phenyl groups were oriented towards the inside of pockets. In fact, the three late states as the next-to-last were in 12 binding pathways, and the total flux of these pathways was up to 52.9 % (26.8 % + 22.6 % + 3.5 % seen in Fig. 3C). In other words, the crystal-like pose of BAY in the end state was directly determined with at least a 52.9 % probability by the inhibitor poses in the three late states (state 22, 21 and 18). Thus, the crystal pose of BAY had at least a 52.9 % probability of being directly converted from poses containing the hydrogen bonds described above.

The hydrogen bonds between E846 and BAY also existed in state 23. Unlike other late states, BAY in state 23 used the metoxyphenyl group to form hydrogen bonds with E846 and M847. The effective flux of this state was 25.2 %, and its only downstream state was the end state (Fig. 3A and 3C). Therefore, the simulated poses of BAY in the end state was 25.2 % likely to be directly converted from the poses in state 23. Considering all cases of late states, the hydrogen bonds between E846 and BAY may determine the simulated poses of BAY in the end state to a 78.1 % extent (26.8 % + 22.6 % + 3.5 % + 25.2 % seen in Fig. 3C).

In all 5 late states, the Top10 common residues were conserved F862, strongly similar residues M847, I866, and L770, weakly similar residues S861, and non-conserved residues L858, M845, H865, and E846 (Fig. 4). Notably, the sequence alignment showed that the residue corresponding to S861 is Glycine in the other 3 PDEs (Fig. S1). As a result, the interactions between S861 and BAY, such as the hydrogen bonds between the side chain of S861 and the metoxyphenyl group of BAY in state 21, were unique to the PDE2A system. In other words, S861 in the PDE2A system should act as a non-conserved residue. As mentioned above, the contribution ratio of non-conserved residues to BAY binding energy was less than 20 % (Fig. 3A). So, non-conserved residues L858, M845 and H865 should contribute a small amount of affinity mainly by virtue of their hydrophobic side chains. Therefore, the non-conserved residues in the late states guided the behavior of BAY through the above non-covalent interactions.

In the end state, state 10, the simulated pose of BAY had two main characteristics of its crystal pose: the π-π stacking between F862 and the pyrazolopyrimidinone group, and the phenyl group entering the hydrophobic subpocket. As a result, the BAY pose in state 10 was the closest to the crystal pose of all states (RMSD = 3.51 Å). Three (I866, L809 and L770) of the Top10 residues came from the hydrophobic subpocket (Fig. 4) [17]. This suggested that hydrophobic pockets may be one of the main sources of PDE2A affinity for BAY. This was also consistent with the hypothesis put forward by previous studies. The above results indicated that state 10 basically restored the crystal binding state of BAY.

In summary, according to the spatial locations of BAY, the eight metastates that comprised the inhibitor binding pathways were classified as the early, middle and late states. In these states, non-conserved residues not only contributed to the BAY binding energy, but also assisted BAY conformational adjustment through non-covalent interactions. These non-conserved residues included R762, I776, R775, I764, Q794, Q781, P869, F777, E846, L858, H865, and so on. Among them, Q794, R762, Q781, and E846 formed hydrogen bonds with BAY, which made BAY form appropriate poses in the corresponding metastates during the binding process, and successfully bind to the catalytic site of PDE2A. In particular, the hydrogen bonds between E846 and BAY directly determined the successful binding of BAY to the catalytic site of PDE2A to a degree of 78.1 %.

Binding pathways of BAY in other three PDE systems

Similarly to the PDE2A system, BAY spontaneously bound to the pockets of PDE4B, PDE5, and PDE10A. However, RMSD values showed that the inhibitor could not form the crystal-like pose in these three pockets (Fig. S2). As mentioned above, the reason for this result may also lie in the binding pathways of BAY in these three systems. Thus, similar to the analysis of the PDE2A system, the conformational analyses would be performed for the metastates comprising the dominant binding pathways.

In the PDE4B system (Fig. S4), the binding pathways of BAY started at the back of the pocket (state 7 and 3), passed through the alpha helix H16 (state 6), the region between the alpha helices H15 and H12 (state 18, 5, and 15) and M−loop (state14,20,19 and 13), and.

terminated at the catalyic pocket (state 10). In all these metastates, the contribution ratios of non-conserved residues to BAY free energy were all less than 40 %, and in the three metastates with RMSD less than 12 Å (state 20,19, and 13), their ratios were less than 20 %. At the same time, except for the two metastates (state 3, 6) far from the pocket, no hydrogen bond were formed between the side chains of non-conserved residues and BAY (Fig. S5). The above data indicated that in the PDE4B system, non-conserved residues did not dominate the affinity for BAY, and no hydrogen bond was formed above M−loop to assist BAY conformational adjustment. As a result, the phenyl group of BAY did not enter the hydrophobic pocket in the end state, and thus this inhibitor did not form the crystal-like pose.

In the PDE5 system (Fig. S6), the binding pathways of BAY started at the region between the alpha helices H15 and H12 (state 18 and 10), passed through the M−loop (state 14,19, and 13), and l-loop (state 16 and 11), and terminated at the catalytic pocket (state 6). In all these metastates, the contribution ratios of non-conserved residues to BAY free energy were all less than 40 %, and in the 5 metastates with RMSD less than 15 Å (state 14,19,13,16 and 6), their ratios decreased to less than 20 %. At the same time, except for the state 11, no hydrogen bond was formed between non-conserved residues and BAY (Fig. S7). The above data indicated that in the PDE5 system, non-conserved residues did not dominate the affinity for BAY, and no hydrogen bond was formed above M−loop to assist BAY conformational adjustment. As a result, BAY did not form the crystal-like pose in the end state; that is, the phenyl group did not enter the hydrophobic subpocket, and the pyrazolopyrimidinone group did not form the π-π stacking with F820.

In the PDE10A system (Fig. S8), the binding pathways of BAY started at the region between the alpha helices H11 and H13 (state 19), passed through the M−loop (state 15, 9, 17, 12, 4 and 16), and terminated at the pocket (state 2). In all these metastates, the contribution ratios of non-conserved residues to BAY free energy were all less than 20 %. At the same time, no hydrogen bond was formed between non-conserved residues and BAY (Fig. S9). The above data indicated that in the PDE10A system, non-conserved residues did not dominate the affinity for BAY, and no hydrogen bond was formed above M−loop to assist BAY conformational adjustment. As a result, BAY did not form the crystal pose in the end state; that is, the phenyl group did not enter the hydrophobic subpocket, and the pyrazolopyrimidinone group did not form the π-π stacking with F729.

In conclusion, BAY moved toward pockets along different binding pathways on the surfaces of PDE4B, PDE5 and PDE10A. In these binding pathways, the non-conserved residues did not dominate the affinity for BAY, and there was no residue that could assist BAY in conformational adjustment like E846 in the system of PDE2A. As a result, BAY failed to adopt a pose that approximates the crystal structure in the catalytic sites of these three PDEs. In other words, the non-conserved residues, especially E846, in the PDE2A system played a important role in assisting BAY conformational adjustment, leading BAY to bind into the catalytic site. This may determine the selective binding of BAY to PDE2A.

Binding pathways of FKG and MC2 in PDE2A system

As mentioned above, the hydrogen bonds between E846 and pyrazolopyrimidinone group may be critical for selectivity of BAY. On the other hand, there are many selective inhibitors targeting PDE2A, and some of them are pyrazolopyrimidinone derivatives [56], [57], [58], [59]. Naturally, an interesting hypothesis is that those hydrogen bonds may also be present in the binding processes of other selective pyrazolopyrimidinone inhibitors targeting PDE2A, such as FKG (IC50 = 2 nM) and MC2 (IC50 = 1.3 nM). Moreover, LiGaMD could be a more efficient way to reproduce the process of inhibitor unbinding and rebinding than unbiased MD due to the potential boost. Therefore, to test the universality of the above hydrogen bonds, LiGaMD was used to reproduce the unbinding and rebinding processes of FKG and MC2 in the PDE2A system. [5].

FKG unbound from the pocket (RMSD greater than 15 Å) and rebound (RMSD less than 5 Å) into the catalytic site of PDE2A 15 times in the about 9.4 μs LiGaMD trajectories (Fig. S3). The minimum RMSD of the simulated poses of FKG reached 1.36 Å. These data indicated that LiGaMD successfully reproduced the unbinding and rebinding processes of FKG in PDE2A system. After the discretization and MSM analysis, the trajectories above were divided into 20 metastates. The state 20 with the lowest RMSD was determined as the end state. Among the four metastates with large RMSDs, state 1 with the smallest flux was set as the initial state. From state 1 to state 20, there were 138 possible binding pathways, among which the flux sum of the first 9 pathways was 60.7 %. In other words, FKG had a 60.7 % chance of successfully binding into the catalytic site along these 9 dominant pathways. Moreover, the metastates making up these dominant pathways participated in 120 of the 138 binding pathways mentioned above, and the total flux of these 120 pathways was 77.8 % (Fig. 5C). That is, FKG had 77.8 % probability of successfully binding to the catalytic site through these metastates. Therefore, the metastates that comprised these dominant pathways should contain information about the molecular determinants of FKG binding selectivity.

Fig. 5

(A) Effective flux (left bar), energy contribution from non-conserved residues (right bar) and RMSD (broken line) of 10 metastates comprising the top 9 binding pathways that dominate the association processes of FKG in PDE2A system. (B) Center of mass (COM) of FKG in those 10 metastates. COMs are displayed as spheres and their colors follow the bars in (A). (C) Flux network from the initial state (state 1) to the end state (state 20). Colors of metastates and their fluxes follow the bars in (A). Non-dominant states are displayed in grey.

There are 10 metastates that made up these pathways (state 1, 3, 11, 7, 19, 15, 8, 17,13, 20, in descending RMSD order). As shown in Fig. 5A, the contribution ratios of non-conserved residues to FKG binding free energy were less than 20 % in all the 10 states. This meant that non-conserved residues did not dominate affinity for FKG.

In addition to affinity, the binding conformational preference of FKG in various states should be another factor affecting its biological properties, which in turn affected its behaviors. In the initial state, FKG bound at the region between alpha helices H15 and H12 (Fig. 5B). This meant that this region was the most likely starting point for FKG binding pathways. Nine of Top10 residues in the initial state were hydrophobic, including two non-conserved residues, H865 and F777. H865 contributed the most to binding energy of FKG (Fig. 6). The hydrophobic interaction of FKG with these residues may result in the inhibitor taking its trifluoromethyl-phenyl group towards the entrance. This pose may facilitate the downward movement of FKG to the entrance in subsequent states.

Fig. 6

Nine metastates comprising the top 9 binding pathways that dominate the association processes of FKG in PDE2A system. For each metastate, upper is the Top10 residues with their contribution to FKG binding energy, and bottom is the poses of FKG (violet stick model) and Top10 residues (stick model) of PDE2A (cartoon model). The conservation of residues is indicated by their colors. Hydrogen bonds are displayed as yellow dotted lines. Metastates are displayed in descending order of RMSD from left to right and from top to bottom. Simulated pose (violet stick model) of BAY is aligned to its crystal pose (cyan stick model) in the end state.

Except for the initial and end states, inhibitors in the other 8 metastates bound near alpha helix H15 at the pocket entrance (Fig. 5B). The total flux of the binding pathways with these 8 metastates is up to 77.8 %. This suggested that these metastates dominated the FKG binding processes. According to the theory of the MSM, the conformation of the end state could only be affected by the conformation of next-to-last state. In fact, the above 8 metastates as the next-to-last state appeared in most binding pathways, and the total flux of these pathsway was up to 76.3 %. In other words, the crystal-like pose in the end state had a 76.3 % probability of being directly converted from the FKG poses in these eight metastates.

Of the above eight states, state 19 had the highest effective flux (Fig. 5A). In this state, S861 on the alpha helix H15 formed hydrogen bonds with the pyrazolopyrimidinone group of FKG (Fig. 6). Although S861 is defined as a weakly similar residue in sequence alignment of 4 PDEs, it is Glycine at this site of the other 3 PDEs (Fig. S1). This determined that the above hydrogen bond could only occur in the PDE2A system. As mentioned earlier, both ligand and receptor need to adopt the appropriate conformation to form hydrogen bonds. Therefore, the hydrogen bond between S861 and pyrazolopyrimidinone groups may play an important role in the conformational adjustment of FKG in state 19.

The hydrogen bond between S861 and FKG also appeared in state 3 (Fig. 6). The total flux of the binding pathways with state 19 and 3 was 30.6 %. These two metastates were the next-to-last states in some pathways whose the flux sum was 24.5 %. According to the MSM theory, these data indicated that FKG had at least a 30.6 % probability of forming a hydrogen bond with S861 by means of its pyrazolopyrimidinone group during its binding, and that the crystal-like pose of FKG had at least a 24.5 % probability of being directly converted from the simulated pose containing the above hydrogen bond.

In state 7, 8, 11, 15, and 17, most of the Top10 residues were conserved or similar residues with hydrophobic sidechains. Among them, the contribution of conserved residue F862 was the largest. In addition, conserved residue F830, strong similar residues M847, I826, I866, L770, weak similar residue S861, and non-conserved residues L858, H865, Y827, T703, M845 also appeared in these states for many times (Fig. 6 and Fig. S10). Therefore, the affinity of FKG on the binding pathways mainly came from the conserved and similar residues represented by F862, but the affinity contributed by the above five non-conserved residues should be also one of the key factors affecting the conformational preference of FKG. As shown in the Fig. 6, in these 5 metastates, inhibitors were located near the alpha helix H15 and took the poses of pyrazolopyrimidinone and trifluoromethyl-phenyl groups pointing out and in pocket of PDE2A, respectively. Among all binding pathways of FKG, the probability of these 5 states being the next-to-last state was 36 %; that is, the probability of the crystal-like pose of FKG being directly converted from the above poses was 36 %.

In contrast to the above poses, the pyrazolopyrimidinone group of FKG in state 13 was inward, while its trichlorogroup was outward (Fig. 6). Compared with the crystal-like pose in the end state, the pyrazolopyrimidinone group of FKG in this state was very close to forming the π-π stacking with F862, and the trifluoromethyl-phenyl group was also close to extending into the hydrophobic subpocket. Compared with the other 8 metastates, the pose of FKG in state 13 was closest to the crystal pose (RMSD = 6.64 Å; Fig. 5A), so it is also closest to the catalytic site consisting mainly of conserved residues. Therefore, the Top10 residues were mainly composed of conserved and similar residues. Moreover, the data showed that the crystal pose of FKG had a 15.8 % probability of being directly converted from the pose in state 13 (Fig. 5C).

Taken together, the conformational preference of FKG in its binding pathways may be caused by the interactions between FKG and S861, L858, H865, Y827, T703, and M845. The hydrogen bond between E846 and pyrazolopyrimidinone group did not form in any of these 10 metastates.

As shown in the Fig. S3, MC2 unbound from the pocket (RMSD greater than 15 Å) and rebinds (RMSD less than 5 Å) into the catalytic site of PDE2A 9 times in the about 7.7 μs LiGaMD trajectories (Fig. S3). The minimum RMSD of the simulated poses of FKG reached 2.15 Å. These data indicated that LiGaMD successfully reproduced the unbinding and rebinding processes of MC2 in PDE2A system.

After discretization and MSM analysis, the above trajectories could be divided into 20 metastates. According to the definitions in Section 2.5, the binding process of MC2 started at state 1 and ended at state 4. Eighteen of the 20 metastates comprised the dominant binding pathways. In these 18 metastates, the contribution ratios of non-conserved residues to the binding free energy of MC2 were less than 20 % (Fig. S11). Most of the non-conserved residues with high contributions were hydrophobic residues, including L858, M845, M848, Y827, H865, etc. Hydrogen bonds only occurred between the pyrazolopyrimidinone group and weakly similar residues S861 (state 6 seen in Fig. S12) and conserved residue Q859 (state 15 and 4 seen in Fig. S13). Thus, the non-conserved residues of PDE2A did not affect the conformational preference of MC2 during the binding process through hydrogen bonds.

To sum up, non-conserved residues did not dominate the affinity of PDE2A for FKG and MC2, and almost all the non-conserved residues with high affinity contributions were hydrophobic residues. Although FKG and MC2 are both pyrazolopyrimidinone PDE2A inhibitors, they exhibited different binding pathways from BAY. Moreover, they hardly formed hydrogen bonds with non-conserved residues during their binding processeses. Therefore, the non-conserved residues did not depend on hydrogen bonding to affect binding processes and binding poses of FKG and MC2 at the catalytic site.

Discussion

In this study, HTMD was performed to reproduce the spontaneous association of BAY, a highly selective inhibitor, with the pockets of four PDE isoforms. Only along the pathways of the PDE2A system can BAY formed the simulated pose approximating the crystal one at the catalytic site. Moreover, LiGaMD was performed to reproduce the unbinding-rebinding processes of FKG and MC2, the other two highly selective PDE2A inhibitors, in the PDE2A system. After the discretization, MSM analysis and free energy calculation, the results showed that (1) non-conserved residues did not dominate the BAY affinity in the four PDE systems; (2) the non-conserved residues outside the pocket may affect the conformational preference of BAY in the binding processes through the hydrogen bonds with BAY; (3) among them, the hydrogen bonds between E846 and BAY may directly determine the successful binding of BAY to the PDE2A catalytic site to a 78.1 % degree; (4) the absence of hydrogen bonds between non-conserved residues and BAY led to the failure of BAY to bind to the catalytic sites of the other three PDEs; (5) the hydrogen bond may not be the key factor affecting the binding pathways of FKG and MC2 by non-conserved residues.

Based on our previous research, MD simulation studies once again showed that different binding pathways led to different binding poses of BAY at the catalytic sites of PDEs, thus determining the selectivity of this inhibitor. In most of the metastates that comprised these binding pathways, the inhibitors took various poses outside the pockets. Meanwhile, the proportion of non-conserved residues outside the catalytic pocket was higher than those inside. Considering the specific biological properties of individual amino acids, non-conserved residues were more likely than other residues to cause differences in binding pathways. Since non-conserved residues did not dominate the affinity of PDE for BAY, this study focused on the conformational preferences of inhibitors determined by the interactions of non-conserved residues with inhibitors. In this way, we found that the hydrogen bonds between E846 and BAY may determine the selectivity of this inhibitor. So, these results suggested that the binding processes of inhibitors, especially the role of non-conserved residues, should be considered more comprehensively in future inhibitor design.

Although BAY, FKG and MC2 are all pyrazolopyrimidinone PDE2A inhibitors, the biggest difference between BAY and the latter two is that BAY has the metoxyphenyl group. As described in Section 3.1, this group formed hydrogen bonds with non-conserved residues Q781 and K872 in state 14, which affected the binding behaviors of BAY in the downstream metastates to a degree of 65.3 %. These binding behaviors included the hydrogen bonds between E846 and BAY, and the conformational preferences that resulted from these hydrogen bonds. The metoxyphenyl group also formed hydrogen bonds with the non-conserved residue E846 and the weakly similar residue S861 in state 23 and 21, respectively. These two metastates directly accounted for 52 % of the successful binding of BAY to the catalytic site. By interacting with the non-conserved residues of PDE2A, the metoxyphenyl group could not only promote hydrogen bonds between E846 and the pyrazolopyrimidinone group, but also directly form the hydrogen bond with E846. According to Section 3.1, the hydrogen bonds between BAY and E846 directly determined the successful binding of BAY at the catalytic site to a total extent of 78.1 %. Therefore, the hydrogen bonds between metoxyphenyl group and non-conserved residues were very important for the successful binding of BAY to the catalytic site of PDE2A, and may also be the reason why the binding pathways of BAY were different from those of FKG and MC2.

Another example of the metoxyphenyl group affecting BAY binding behavior comes from LiGaMD of PDE2A-BAY system. The RMSD of BAY was consistently lower than 15 Å in the 7.6 μs LiGaMD (Fig. S3). In other words, BAY failed to run out of the PDE2A pocket even with the addition of the potential boost. However, FKG and MC2 lacking the metoxyphenyl group could escape and enter the pocket of PDE2A many times at the same time scale. Therefore, BAY should have a longer residence time than FKG and MC2. As shown in Fig S14, at least eight of the 14 metastates that made up the dominant pathways showed the metoxyphenyl group oriented toward the inner side of the PDE2A pocket dominated by hydrophobic residues. Therefore, it was speculated that the interactions between the metoxyphenyl and the hydrophobic residues inside the pocket should retards the release BAY from from the PDE2A pocket.

In this study, the most important basis to determine the key interactions were the conversion rates between the metastates, namely fluxes. Although sophisticated cluster analysis can detect metastates [60], the fluxes relied on the MSM analysis. By means of the MSM, the initial and end states of each system were determined, and the analysis scope was limited to the binding pathways whose the sum of fluxes was just more than 60 %. In PDE-inhitiors systems, the metastates that comprised these binding pathways were involved in most of inhibitor binding events (flux sum greater than 77.8 %). This validated the above limits of the analysis scope and implied that other metastates may have little to do with the binding processes of inhibitors. Further flux and conformational analysis identified key interactions that may dominate the inhibitor binding processes and suggested that the selectivity of PDE2A inhibitors may depend on their conformational preferences. These results showed that the combination of the MD simulation and the MSM can play an important role in analyzing the molecular mechanism of ligand-receptor interactions over long-time scales. In fact, the combination of these two methods has been successfully used to explore the molecular determinants of the inhibitor residence time and to analyze the allosteric effects between antagonists [61], [62].

In conclusion, the non-covalent interactions between BAY and non-conserved residues determines the conformational preference of the inhibitor in its binding process, and thus determines the binding pathways of BAY on the surface of PDE isoforms, which is ultimately manifested as the high selectivity of BAY to PDE2A. Among them, the non-covalent interactions of the metoxyphenyl group with the non-conserved residue of PDE2A make the binding processes of BAY different from those of FKG and MC2, and also facilitates the formation of hydrogen bonds between E846 and BAY. These hydrogen bonds directly determine the successful binding of BAY at the catalytic site of PDE2A to a total extent of 78.1 %. Thus, this study demonstrates a possible molecular mechanism by which the binding process affects BAY selectivity, and thus facilitates the further development of BAY-derived inhibitors targeting PDE2A.

CRediT authorship contribution statement

Qing Liu: Investigation, Formal analysis, Visualization, Writing – original draft, Writing – review & editing, Conceptualization. Menghua Song: Investigation, Writing – review & editing. Yue Qiu: Investigation, Writing – review & editing. Elaine Lai-Han Leung: Supervision, Conceptualization, Writing – review & editing. Qiang Huang: Supervision, Conceptualization, Writing – review & editing. Xiaojun Yao: Supervision, Conceptualization, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.