Received Date: October 30, 2015; Accepted Date: November 02, 2015; Published Date: November 05, 2015
Citation: Natarajan P, Swargam S, Hema K, Vengamma B, Umamaheswari A (2015) E-Pharmacophore Based Virtual Screening to Identify Agonist for PKA-Cα. Biochem Anal Biochem 4:222. doi:10.4172/2161-1009.1000222
Copyright: © 2015 Natarajan P, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Biochemistry & Analytical Biochemistry
Owing to PKA-Cα unique functions in regulating tau splicing alternatively in neurons results in aggregation of tau which contributes to neurofibrillary tangles and taupathies generation a hallmark of Alzheimer’s disease (AD). PKACα also inhibits GSK3β thus it has been a therapeutic target for AD intervention. In this study, e-pharmacophore and multiple docking strategies were followed to propose a novel PKA-Cα agonists. Nine e-pharmacophores were developed from nine co-crystal structures such that all the critical pharmacophoric features involved in their bioactivity of PKA-Cα were effectively mapped. Rigid receptor docking (RRD) was performed with the library of PKA-Cα activators having 3512 shape screened compounds towards PKA-Cα. To derive the best leads, dock complexes were further subjected to QPLD, IFD and MM-GBSA calculations. PKA-Cα-lead1 dock complex was subjected to 50 ns MD simulations run. Comparative analysis between obtained 25 leads and 9 co-crystal ligands revealed three the best leads. Among the three, lead1 has the least docking score with lowest binding free energy with better binding orientation towards PKA-Cα. The constancy of PKA-Cα-lead1 interactions was revealed by 50 ns MD simulations run. Thus ADME predictions and results from RRD, QPLD, IFD and MD simulations affirmed that the proposed three leads might be used as a potent agonists for PKA-Cα.
PKA-Cα; 3R-tau; 4R-tau; e-pharmacophore; Alzheimer’s disease; MD simulations
Alzheimer’s disease (AD) is a common neuro degenerative disease characterized by a slow, progressive decline in cognitive function and behavior. Along with the deposition of Amyloid plaques, formation of neurofibrillary tangles and tauopathies are the hallmarks of AD . Tau is defined with stabilizing the microtubules assembly which gives definite structure to the cell. For normal functioning of brain the ratio of 3R-tau to 4R-tau should be balanced that is to be 2:1 ratio, but cAMPdependent protein kinase catalytic subunit alpha (PKA-Cα) deregulates the tau splicing alternatively .
The down regulated of PKA-Cα deregulates tau splicing alternatively which results in elevated levels of 3R-tau production by excluding the exon 10. In the increased 3R tau levels microtubules dissociates with the protein tau, thus hydrophobic tau tends to oligomerization resulting in aggregation of tau oligomers. This ensues the hyper phosphorylation of tau oligomers, which contributes neurofibrillary tangles in the hippocampus region of the brain leading to AD . Lace et al., also reported that down regulation of PKA-Cα correlates with the increase in 3R-tau expression levels . Down regulation of PKA in AD brain correlates with increase in 3R-tau expression. PKA-Cα also inhibits GSK-3β, thus it has been an attractive therapeutic target for AD intervention . Thus alteration in the 3R-tau/4R-tau ratio is sufficient to trigger neurofibrillary tangles formation and loss of neuronal cell structure resulting in neurodegeneration of hippocampus region of brain in frontotemporal dementia.
Dysregulation of exon 10 causes neurofibrillary degeneration by destabilizing the microtubules associated with tau. PKA-Cα has emerged as a key kinase that is able to interact with many of the proteins involved in the etiology of AD as well as other tauopathies . Thus by elevating the PKA-Cα levels may result in reducing the 3R-tau/4R-tau ratio and also PKA-Cα inhibits GSK-3β. As there is a wealth of structural information of PKA-Cα available, e-pharmacophore mapping and its use in e-pharmacophore-based virtual screening and multiple docking strategies were carried out in the present study.
In the current workflow, all the analysis were executed in HPZ800 workstation (12 processors; 12 GB RAM and Sun Microsystem workstation running on Cent OS 5. Receptor-based pharmacophore approach was performed due to the availability of twenty two crystal structures of the PKA-Cα. Co-crystalized structures having resolution ≤ 2.0 Å (nine) were considered for the study and which affords a range for search of novel scaffold or compound that could elevate the binding site of PKA-Cα.
Active site analysis and rotein structure preparation
Nine crystal structures such as 3OVV, 3POO, 3AMA, 3OWP, 3L9M, 3VQH, 3OOG, 3L9L, 3L9N with a resolution less than 2.0 Å crystalized with the co-crystal ligands 1SB, S69, SKE, 2SB, L9M, IQB, YTP, L9L, L9N respectively were retrieved from the protein databank (PDB) . Active site was predefined for further study as these residues contribute to the structural and functional properties of the protein. The co-crystal ligands (substrates or activator) interactions towards PKA-Cα were analyzed using PDBsum .
Using the OPLS_2005 force field, by converging the heavy atoms to RMSD of 0.3 Å protein structure minimization was carried out. Hydrogen atoms were added to all the atoms in the system and bond orders, formal charges were added for the hetero groups. Optimization was done at neutral pH [7-9].
Generation of e-pharmacophore hypothesis
Nine co-crystal ligands were prepared using LigPrep v.3.0 (Schrodinger) and were docked into their respective crystal structures using Glide XP docking. Nine refined complexes were selected for e-pharmacophore modeling. The refined ligand-receptor complexes along with Glide XP descriptor information, the energy terms were mapped on to the atoms [7,10,11]. Glide XP energies were added together for all the atoms that comprises the pharmacophoric site. As fitness score is the measure of how well the ligand fits into the receptor with the reference to the ligand interaction energies, pharmacophoric sites with fitness score less than -0.5 were rejected [7,10].
Pharmacophore model evaluation
Nine e-pharmacophore models of PKA-Cα were generated from nine co-crystal ligands. The nine e-pharmacophore hypotheses were considered for validating the virtual screening protocol by using enrichment factor (EF) metrics calculations. The co-crystal ligands were considered as actives which were then combined to decoys set having 1000 drug-like compounds retrieved from Schrodinger  to form an internal library of 1009 compounds. The ability of the e-pharmacophores to differentiate the actives from internal library was evaluated for validating the screening protocol. The nine e-pharmacophores were taken as a query to screen the internal library to generate a library of PKA-Cα activators and was docked to PKA-Cα. The results availed were analyzed with enrichment factor at 1%, which measures the enrichment for recovering the known actives and Boltzmann-enhanced discrimination of receiver operating characteristic curve (ROC) i.e., BEDROC (α=20) metrics which measures the enrichment for early recognition of actives among the ranked internal library .
E-pharmacophore based database screening
E-pharmacophores that are validated were considered for shape screening. Shape based similarity screening was performed using PHASE module of Schrodinger for the multiple e-pharmacophore models generated with the selected features [7-9]. Flexible search was carried out during the shape based similarity screening against more than one million compounds derived from small molecule databases such as ChemBank, ChemPDB, KEGG ligand, Anti-HIV NCI, Druglikeness NCI, Not annotated NCI, AKos GmbH and Asinex Ltd. Distance matching tolerance of 3.0 Å were considered for the hits matched for the pharmacophoric features [14-17]. Multiple conformers were generated for all the tautomer’s matched from the screening. In order to reduce the false positives, inactive compounds were rejected. To refine the screening process receptor-based excluded volumes were also integrated in shape screening. All the hits obtained from the shape based similarity search based on the e-pharmacophore models and nine co-crystal ligands were exported as library of PKA-Cα activators .
Three tier docking such as rigid receptor docking (RRD), quantum polarized ligand docking (QPLD) and induced fit docking (IFD) protocols were employed to predict the scoring and binding interactions between PKA-Cα and the ligands. The prepared in-house library was then docked into the PKA-Cα active site. The co-crystal ligands taken for the development of e-pharmacophore models were docked with the same protocol as it can be a basis for comparison [7-9].
Using Epik  and LigPrep all the ligands were prepared by accomplishing the chemical correctness with protonation, stereo chemical and ionization variations. Energy minimization was done at neutral pH 7.0 ± 2.0 units . Ligands with reactive functional group, high ionization energy/tautomer states were removed from the generated conformations. Ligands that are not obeying the Lipinski’s rule of five were discarded from the multiple conformations generated . A receptor grid of 10x10x10 Å was generated around the active site residues of PKA-Cα (3OVV) crystal structure using Glide v.5.9 (Grid based Ligand Docking with Energetics).
The PKA-Cα grid and library of PKA-Cα activator compounds were exported to Maestro v.9.6. . RRD protocol has a three level docking such as high throughput virtual screening (HTVS), standard precision (SP) and extra precision (XP) to screen the ligands selectively from lesser stringency to higher stringency at every level of glide docking [7,8,22]. Finally the XP docked complexes were evaluated using XP Glide Score (XPG Score) and best pose was generated as output, further taken for free binding energy (ΔG) calculations [7-9]. XPG Score is an empherical scoring function that approximates the ligand binding energy and the parameters such as force fields, penalties for the interactions that have the influence of ligand binding with the receptor. The XPG Score is given by
GScore = a*vdW+ b*Coul + Lipo + Hbond + Metal + BuryP + RotB + Site
Whereas vdW denotes van der Waals energy, Coul denotes Coulomb energy, Lipo denotes lipophilic contacts, HBond indicates hydrogen-bonding, Metal indicates metal-binding, BuryP indicates penalty for buried polar groups, RotB indicates penalty for freezing the rotatable bonds, Site denotes polar interactions with the residues in the active site and the a = 0.065 and b=0.130 are coefficient constants of van der Waals energy and Coulomb energy respectively.
Quantum polarized ligand docking (QPLD) was employed to compute the atomic partial charges of the leads through quantum mechanical and molecular mechanical (QM/MM) calculations . Q-Site module of Schrodinger was employed for QPLD docking . Initial docking was carried on using Glide standard precision mode to generate five diverse poses for each lead. Atomic partial charges for PKA-Cα-lead complexes were calculated based on the proposed orientation of leads relative to the PKA-Cα and influenced by the electrostatic field of the PKA-Cα. Single point energy determination was treated by the QPLD to determine the most energetically favorable lead poses with respect to PKA-Cα. Intermediate charge sets of PKACα and lead complexes calculated earlier in the Q-Site refinement were re-docked through Glide extra precision mode to get the empirical atomic charges and further carried out for binding free energy (ΔG) calculations .
Active site residues and their flexibility of PKA-Cα were considered for IFD protocol in Schrodinger. Leads on interacting with binding site residues of PKA-Cα, it undergoes side-chain or backbone conformational change or both. These conformational changes allow the PKA-Cα to generate closest conformer to the shape and binding mode of the leads. IFD was carried out with default parameters such as VdW radii was scaled to 0.5 Å and 20 conformational poses were calculated for both PKA-Cα and lead atoms . OPLS_AA 2005 force field with implicit solvation model was applied for energy minimization process. Maestro (Schrodinger) interface and Python scripts automatize the IFD method in three consecutive steps. The active site and dimensions for the energy grid was defined around co-crystal ligand for initial docking using Glide XP mode .
Later by using Prime v.3.6, PKA-Cα-lead complexes side chain and backbone refinement was performed to each docked conformations. The PKA-Cα-lead complexes were ranked by Prime energy. Receptor structures below 30 kcal/mol of the minimum energy structure were passed through for a final round of Glide docking and scoring. In the last phase, Glide XP was utilized to re-dock each leads into every refined low-energy conformation of PKA-Cα produced in the earlier phase. Side chain orientations were automatized in IFD with the insertion of Prime . PKA-Cα-lead interaction energies and total energy of the system was calculated as IFD score. Based on the IFD scores the poses generated were ranked [2,7,26].
Free energy calculations
Using molecular mechanics/generalized Born surface area (MM/ GBSA) the binding free energy (ΔG) of PKA-Cα-lead complexes were calculated by Prime approach . For each PKA-Cα-lead complexes ΔG was calculated by using the equation as follows:
ΔGbinding = ΔGcomplex – (ΔG PKA-Cα + ΔGlead)
WhereΔGcomplex, (ΔG PKA-Cαand ΔGlead are the free energies of complex, PKA-Cα and leads respectively.
Through MM/GBSA OPLS_2005 molecular energies (EMM), polar solvation through surface generalized Born solvation model (GSGB) and a nonpolar solvation term (GNP) composed of nonpolar solvent accessible and vdW interactions were calculated. The binding free energy (ΔG) calculations are much more accurate than the XPG Score [28,29]. ΔG calculations can be given more elaborately as follows
ΔGbind = ΔE + ΔGsolv + ΔGSA
Where, ΔE is the minimized energies, ΔGsolv solvation free energies, ΔGSA is the difference in surface area energy of the PKA-Cα- lead complex and sum of the surface energies of PKA-Cα and leads respectively.
Simulations were carried out through Prime v.3.6 by implementing generalized Born surface area (GBSA) continuum model . Linear combination of pairwise overlap method was used to determine the solvent-accessible surface area. Instead of vdW surface Gaussian surface was employed for better representation of the solvent accessible surface area [28,29].
Molecular dynamics simulations
Inter-molecular and intra-molecular interactions that contributes the complex stability can be analyzed using molecular dynamics (MD) simulations. To evaluate the stability of the PKA-Cα-lead1 dock complex molecular dynamics calculations [7-9,29] were executed for 50 ns through Desmond v3.3. The PKA-Cα-lead1 dock complex was solvated using orthorhombic simple point charge (SPC) water model. Counter ions were used to neutralize the solvated system and physiological salt concentration was limited to 0.15 M. OPLS-AA force field was consigned to the solvated dock complex system [30,31]. Periodic boundary conditions were specified to the solvated system. Lennard-Jones interactions cut off were limited to 10 Å and particle mesh Ewald (PME)  method was employed to treat the electrostatic interactions. SHAKE algorithm  was implemented for limiting the movements of all the hydrogen atoms involved in the covalent bonding.
Multistep protocol was employed to simulate the PKA-Cα-lead1 dock complex solvated system in Maestro v.9.6 proceeding to the molecular dynamics simulations production run. In the earlier stage only solvent molecules were allowed to minimize, while the PKA-Cα and lead1 are kept fixed. Then the entire system was minimized in the later stages using hybrid method of steepest descent and limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm [21,33-35]. Maximum of 2000 steps with solute restrains and convergence criteria was set to 50 kcal mol-1Å-1 followed by energy minimization for 2000 steps without solute restraints and the convergence criteria was limited to 5 kcal mol-1Å-1. For restraining the non-hydrogen solute atoms, at a temperature of 10 K with a thermostat relaxation constant of 0.1 ps a short 12 ps simulation was carried out in the NVT ensemble. Followed by short 12 ps simulation in NPT ensemble using 10 K temperature (with thermostat relaxation constant of 0.1 ps and barostat relaxation constant of 50 ps) for restraining the non-hydrogen solute atoms. Then in the NPT ensemble solute non-hydrogen atoms were restrained for 24 ps simulations at a temperature of 300 K (with a thermostat relaxation constant of 0.1 ps; barostat relaxation constant of 50 ps) and 24 ps simulations in the NPT ensemble were carried out with no restrains at 300 K temperature (with a thermostat relaxation constant of 0.1 ps; barostat relaxation constant of 2.0 ps). The solvated system temperature and pressures were controlled by using Berendsen thermostat and barostat [33-35].
A multiple time step RESPA integration algorithm was used for bonded, near non-bonded and far non-bonded interactions in the dynamics production run. Following the relaxation of solvated system, a 50 ns of dynamics production run in the NPT ensemble (at temperature of 300 K, with thermostat relaxation time of 1.0 ps; 1.01 bar pressure, with barostat relaxation time of 2.0 ps) using a Nose-Hoover thermostat and Martyna-Tobias-Klein barostat was performed [36,37]. Energy and trajectory atomic coordinate data were recorded with a time interval of 4.8 ps. The monotony of PKA-Cα-lead1 dock complex was sound out by energy potential, root mean square deviations (RMSD), root mean square fluctuations (RMSF), inter-molecular hydrogen bond interactions monitoring.
PKA-Cα belongs to the class of kinase family. Among the twenty two co-crystal structures of PKA-Cα, nine co-crystal structures that are having resolution of ≤ 2.0 Å in the PDB were considered for the present study.
Protein structure preparation and active site analysis:
The resolved crystal structures of PKA-Cα were reported with nine β-sheets and fifteen core helices. The nine co-crystal structures with ≤ 2.0 Å resolution were retrieved from the PDB. The active site defines the favorable surface which were desirable for ligand binding towards receptor hence the active site was analyzed. The active site was defined with the residues present around 4 Å region surrounding the ligand which was co-crystalized in PDB structures. Pre-defined active site was cross checked with PDBsum. As PKA-Cα is a kinase, the active site residues resides deep in the cleft between the β-sheets and helices along with the allosteric site residues that are responsible for PKA-Cα binding with the ligands.
As multiple resolved structures were available for PKA-Cα in the PDB, energy-optimized structure based pharmacophore method was practiced in the present study. E-Pharmacophores were developed for nine co-crystal structures of PKA-Cα by Glide XP docking. As pharmacophoric site based on the structural and interactional energy information between the PKA-Cα and the co-crystal ligands, nine e-pharmacophores were developed from nine co-crystal ligands such that all the Glide XP energetic terms were mapped on to the atoms. E-Pharmacophores were written with the selected features such that, it could effectively map all the pharmacophoric features which were responsible for PKA-Cα bioactivity.
The generated nine e-pharmacophores were validated using enrichment studies and were used for screening the small molecule databases. The generated nine e-pharmacophores and number of geometrically similar compounds matching to the pharmacophoric features screened were given in the Table 1. Each e-pharmacophores differ from another in screening performance based on the pharmacophoric features. Overall small molecule databases screening efficiency was improved by multiple e-pharmacophores. The derived e-pharmacophores had five, four and three featured models with pharmacophoric sites, such as acceptors (A), donors (D), aromatic ring (R) and hydrophobic group (H).
|S.No.||PDB ID||Resolution in Å||Pharmacophore feature||PKA-Cα e-pharmacophores||No. of Hit compounds|
Table 1: Generated e-pharmacophores of PKA-Cα
All nine PKA-Cα e-pharmacophore hypotheses were generated with five, four and three featured models respectively. SKE co-crystalized with 3AMA, had five featured model with the pharmacophoric features such as ADRDR. Co-crystal ligands 2SB, IQB and L9L of 3OWP, 3VQH and 3L9L respectively were written as four featured models and the pharmacophoric features were given as RRHD, RARR, and ARRH respectively. Remaining five co-crystal ligands of L9M (3L9M), L9N (3L9N), YTP (3OOG), 1SB (3OVV) and S69 (3POO) had three featured models with RRR, RRR, RAH, RAR and RAR pharmacophoric features respectively (Table 1). Features having less than -0.5 fitness score were rejected, as the fitness core defines how well the site point interacts to the receptor in reference structures. The nine different e-pharmacophores developed from nine co-crystal ligand structures of PKA-Cα increases the probability of identifying the diverse and active compounds complementary with the flexibility of residues interaction in the PKA-Cα active site [7,10,11].
Pharmacophore model validation
Known nine co-crystal ligands of PKA-Cα were considered as actives and decoys set having 1000 drug-like compounds from Schrödinger  were combined to form an internal library of 1009 compounds. All the nine e-pharmacophores were scrutinized for the ability of differentiating the actives and decoys from the developed internal library. Using the find match option in shape screening, nine e-pharmacophores were taken as a query to screen against the internal database developed. The 173 hits obtained were considered as PKA-Cα activators internal library, which has retained all the known actives. The generated PKA-Cα activators internal library was then docked into the active site region of PKA-Cα. Commonly used metrics, such as EF1% and BEDROC (α=20) were considered to evaluate the adopted virtual screening protocol .
The enrichment curve graphically represents the quality of known actives that are retrieved from the internal library, which were ranked by comparing to the decoys (Figure 1). EF1% measures the enrichment for recovering the known actives relative to the randomly distributed decoys within the defined internal library. 44% of actives are retrieved within the 1% of enrichment factor from the internal library that comprises both actives and decoys (EF1% = 44). Thus all the 9 actives were retrieved from the 44% of retrieved actives. ROC metric corresponds to the position of actives to the orderly ranked compounds that are linearly arranged among the internal library defined. Truchon and Bayly considered ROC with ≥ 0.7 as a suitable execution measuring value (where, ROC was limited to 0-1) . In the present study all the known actives were retrieved with ROC of 0.99 relative to the PKA-Cα e-pharmacophores in the virtual screening.
As suggested by Truchon and Bayly, BEDROC metric was set to α=20.0 which entails that 80% of the maximum contribution to the BEDROC comes from the first 8% of the ordered list . BEDROC (α=20) metrics measures the early recognition enrichment of actives among the ranked compounds from the internal database. BEDROC value of 0.88 is a beneficial value that embodies the magnitude of early recognition of actives from the ranked compounds in the internal library. Deduced EF, ROC and BEDROC values were 44 (EF1%), 0.99, 0.88 (α =20.0) respectively, elucidates that the nine e-pharmacophores considered for virtual screening against more than one million small molecule were efficacious and sufficient in retrieving the active compounds.
Structurally similar compounds exhibits similar activity, hence structurally and geometrically similar compounds were retrieved form shape screening against small molecule databases having more than one million compounds that are matching to the e-pharmacophores were generated. PKA-Cα activators library having 3503 compounds was developed from the shape screened compounds along with the nine co-crystal ligands. 2888 compounds have passed the Lipinski’s filter. Subsequently 2580 compounds resulted from reactive filters were considered along with the nine co-crystal ligands for further docking studies and free energy calculations.
To study the basis of molecular interactions between PKA-Cα and ligand compounds from the library of PKA-Cα activators and to analyze the binding affinity, all the ligands were docked in to the receptor. To define representative structure among the nine co-crystal structures, RMSD calculations were carried out with the docked complexes. Root mean square deviation between nine co-crystal structures to their docked structures with their respective ligands were computed (Table 2). Among the nine co-crystal structures, 3OVV has the lowest RMSD of 0.00 Å and the lowest crystal structure resolution of 1.58 Å. Thus, the lowest RMSD with lower resolution made 3OVV PDB structure as a representative structure to perform docking studies so as to get a correlation between biological activity and binding free energy towards the different ligands. 2580 ligands docked into the PKA-Cα (3OVV) grid of 10x10x10 Å generated around the centroid of active site resulted 258 compounds from HTVS and 26 compounds from SP which were further submitted to XP which disclosed 25 ligand molecules and the ligand molecules were reassessed for free energy calculations using MM-GBSA. As multiple conformers were generated, to deduce the best scoring leads further they were subjected to QPLD and IFD.
|S.No.||PDB ID||RRD RMSD (Å)||QPLD RMSD (Å)||IFD RMSD (Å)|
Table 2: RMSD of co-crystal ligands with respect to their docked orientation after RRD, QPLD and IFD.
Among the three leads, the lead1 (5-(2-(3-methoxy phenyl) amino) benzimidazol (1,2-a) quinolone -6-carbonitrile) showed the highest binding affinity with ΔG value of -87.77 kcal/mol and XPG Score of -6.92 kcal/mol. Lead1 formed two hydrogen bonds with PKA-Cα. Backbone atom of Phe-54 formed non-covalent hydrogen bond with the oxygen atom of anisole moiety of lead1. Oxygen atom of Lys-72 showed π-cation interaction with the aromatic ring of anisole moiety of lead1. Negatively charged residues such as Glu-121, Glu-127, Asp- 184 and hydrophobic residues such as Leu-49, Val-57, Ala-70, Ile-73, Leu-74, Val-104, Met-120, Tyr-122, Val123, Leu-173 and Phe-327 were involved in van der Waal’s interactions towards lead1. The docking interactions of PKA-Cα-lead1 were depicted in Figure 2.
Pradeep et al., 2015, reported that quantum polarized ligand docking interactions are reliable for ranking the leads than the regular Glide XP scoring . QM was utilized for calculating the atomic charges with OPLS force field by applying semi-empherical method in QPLD. From the analysis of QPLD results for the nine co-crystal ligands disclosed that S69 crystallized with 3POO had the least ΔG score of -78.93 kcal/mol and QPLD score of -6.82 kcal/mol. Backbone atom of Phe-54 formed hydrogen bond with the oxygen atom of anisole moiety of S69. Side chain atoms of Lys-72 formed hydrogen bond with the oxygen atom of aceto-hydrazide moiety of lead1. Backbone atoms of Glu-121 and Val-123 formed hydrogen bonds with S69 (Figure 3a). On comparing the binding free energies of obtained ligands to the co-crystal ligand complexes, lead1 had least binding free energy (ΔG as -85.9 kcal/mol and QPLD score as -6.81 kcal/mol) with similar binding pattern observed in Glide XP docking and found to be stable. Backbone atoms of Phe-54 formed hydrogen bond with the oxygen atom of anisole moiety of lead1. Lys-72 showed π-cation interaction with the aromatic ring of anisole moiety of lead1. Backbone atoms of negatively charged Asp-184 formed hydrogen bond with nitrogen atom of carbonitrile moiety of lead1. Hydrophobic residues such as Leu-49, Val-57, Ala-70, Ile-73, Leu-74, Val-104, Met-120, Tyr-122, Val123 and Negatively charged residues such as Glu-121, Glu-127 along with the positively charged residues such as Arg-56 and Lys-52 were involved in van der Waal’s interactions which stabilizes the PKA-Cα-lead1 complex (Figure 3b). Similar binding interactions were observed in lead docking complexes with the leads docked around the PKA-Cα active site and allosteric site residues.
PKA-Cα showed three hydrogen bonds with the lead1 in IFD docked complex. Lead1 showed IFD score as -6.92 kcal/mol and ΔG value as -87.94 kcal/mol. Lead1 interactions in IFD docking mode with PKA-Cα are shown in Figure 4. Backbone atoms of Phe-54 formed hydrogen bond with the oxygen atom of anisole moiety of lead1. Side chain atoms of Lys-72 formed hydrogen bond with nitrogen atom of carbonitrile moiety of lead1. Lys-72 also exhibited π-cation interaction with the aromatic ring of anisole moiety of lead1.
Pradeep et al, reported that absolute closer experimental affinities can be derived by the including the entropic effect on the binding free energy (ΔG) calculated using MM-GBSA . The binding free energy (ΔG) derived through MM-GBSA calculations had much more consistency than the docking scores and improves the ranking of potential leads as suggested by Yang et al., . Similar approach was practiced in the present study and similar results were obtained in XP, QPLD and IFD. Three tier docking strategy resulted three best leads for activating the PKA-Cα. Different principal descriptors and ADME properties of three proposed agonists and co-crystal ligands were calculated and tabulated in Table 3 and 4.
|Lead||MW||Rotor||SASA||FOSA||WPSA||PISA||Volume||Donor HB||Accept HB||IP||EA||Glob|
Foot note: (Range 95% of Drugs)
(Range 95% of Drugs) (130.0 / 725.0)
MW = Molecular Weight (130.0/725.0) (0.0/15.0)
Rotor=No. of Rotatable Bonds (0.0/15.0) (1.0 / 12.5)
Dipole=Dipole Moment (1.0/12.5) (300.0 /1000.0)
SASA =Total solvent accessible surface area (300.0/1000.0) (0.0 / 750.0)
FOSA =Hydrophobic solvent accessible surface area (0.0/750.0) (7.0 / 330.0)
FISA =Hydrophilic solvent accessible surface area (7.0 /330.0) (0.0 / 450.0)
PISA =Carbon Pi solvent accessible surface area (0.0/450.0) (0.0 / 175.0)
WPSA=Weakly Polar solvent accessible surface area (0.0/175.0) (7.0 / 200.0)
PSA=vdW Polar surface area (7.0/200.0) (500.0 /2000.0)
Volume=Molecular Volume (A^3) (500.0/2000.0) (0.0 / 6.0)
Donor= Donor - Hydrogen Bonds (0.0/6.0) (2.0 / 20.0)
Accpt HB= Acceptor - Hydrogen Bonds (2.0 /20.0)
IP (eV) =Ionization Potential (7.9 /10.5) (7.9 / 10.5)
EA (eV) =Electron Affinity (-0.9 /1.7) (-0.9 / 1.7)
Glob=Globularity (0.75 /0.95) (0.75 / 0.95)
The abbreviations used to specify different principal descriptors of three proposed agonists in the table and their range in 95% of the available drugs are given. The range for properties of 95% drug is given based on QikProp.
Table 3: Principle descriptors of proposed PKA-Cα agonists and co-crystal ligands.
|Leads||Log Po/w||Log S||CIlog S||Log BB||Rule of 5||Rule of 3||Log Kp||Log KhSa|
Foot note (Range 95% of Drugs)
LogP o/w = log P for octavo/water (-2.0/ 6.5) (-2.0/ 6.5)
Logs = log S for aqueous solubility (-6.5/0.5) (-6.5/0.5)
CIlogS = log S - conformation independent (-6.5/0.5) (-6.5/0.5)
LogBB = log BB for brain/blood (-3.0/1.2) (-3.0/1.2)
Log KP = log KP for skin permeability (KP in cm/hr) (KP in cm/hr)
Log Khsa = log K hsa Serum Protein Binding (-2.5/1.5) (-2.5/1.5)
Lipinski Rule of 5 Violations (maximum is 4) (maximum is 4)
Jorgensen Rule of 3 Violations (maximum is 3)
Table 4: ADME properties of proposed PKA-Cα agonists and co-crystal ligands.
Molecular dynamics simulations
MD simulations were performed to know the inter-atomic interactions that facilitate the complex stability. To analyze the conformational stability and steady nature of PKA-Cα-lead1 complex in the physiological environmental conditions embedded with water molecules, temperature and pressure molecular dynamics simulations studies were performed up to 50 ns (nanoseconds). For comparative analysis, co-crystal ligand S69 complex was also simulated for 50 ns. Lys-72 had shown ~20% of π-cation interaction and backbone atoms had shown ~25% of hydrogen bonding interaction with S69 during the 50 ns simulations run. Backbone atoms of Glu-121 formed 78% of hydrogen bond with the hydroxyl atom of phenol moiety of S69 in all the 10416 trajectories. Backbone atoms of Val-123 formed ~48% of hydrogen bond with the hydroxyl group of phenol moiety of S69.
Analysis of PKA-Cα-lead1 docked complex MD simulations revealed a stable interaction pattern and bonding network than the co-crystal ligand S69. The energy plot exhibited the consistency of the PKA-Cα-lead1 docking complex during 50 ns MD simulations run in the acceptable solvated model condition (Figure 5). RMSD was calculated for defining the resulted structural rearrangements in forming a PKA-Cα-lead1 stable complex during 50 ns simulations run. An acceptable range of root mean square deviation (RMSD) of PKA-Cα-lead1 complex is observed from 10416 trajectories (Figure 6). RMSD range for PKA-Cα backbone and lead1 were 0.80 Å to 1.80 Å and 0.60 Å to 1.00 Å respectively and found to be in the conservative range, which is stable throughout all trajectories. Average RMSD for PKA-Cα backbone was 1.20 Å and an average RMSD range for lead1 was 0.50 Å. Root mean square fluctuation (RMSF) was used to calculate the fluctuations of backbone and side chain residues of PKA-Cα during 50 ns simulation period. The RMSF was calculated by averaging all atoms of the given residues in the MD trajectories. RMSF for binding site residues of PKA-Cα were within the limit of 4.0 Å. Average RMSF for PKA-Cα backbone and side chain were 1.20 Å and 1.60 Å respectively (Figure 7). The lower RMSD of PKA-Cα-lead1 complex shows that lesser rearrangement of PKA-Cα towards lead1. The lower RMSF value of the PKA-Cα-lead1 complex indicates the reduced random motions and minimal fluctuations of the backbone residues to the side chain residues during the 50 ns simulations run. Therefore, lower RMSD and RMSF values of PKA-Cα-lead1 complex indicate the lesser structural re-arrangements and reduced internal motions in and around the active site and allosteric site residues, respectively.
Backbone atoms of Phe-54 formed hydrogen bond with the oxygen atom of anisole moiety of lead1. Phe-54 also showed π- π interaction with aromatic ring of anisole moiety of lead1. Lys-72 showed π-cation interaction with the anisole moiety of lead1. Lys-72 formed salt bridge with the nitrogen atom of carbonitrile moiety of lead1. Hydrophobic residues such as Leu-49, Val-57, Ala-70, Leu-74, Tyr-122, Val-123, Leu-173 and Phe-327 were involved in Van der Waal’s interaction towards PKA-Cα. During 50 ns simulations run lead1 formed similar interaction pattern towards PKA-Cα as observed in QPLD and RRD. Water molecules also facilitated the bonding of PKA-Cα with lead1 in the 50 ns simulations run (Figure 8). Additionally, water bridges were formed with Glu170, Asn-171 and Asp-184 of PKA-Cα. Poornima et al, reported that electrostatic interactions can be increased by 10 kcal/mol with the increase of receptor-ligand water bridges . Thus additional displacement of ordered water molecules increases the binding affinity and entropy gain in the complex formation [40-42]. Bonding interactions of active site residues were observed to be stable in all the 10416 trajectories and the complexes were found to be consistent. Phe- 54 formed ~46% of hydrogen bond with oxygen atom of anisole moiety of lead1 and ~75% with aromatic ring of anisole moiety of lead1 in all the 10416 trajectories. Lys-72 formed hydrogen bond with lead1 in ~76% with anisole moiety and ~76% with carbonitrile moiety of lead1 in all the 10416 trajectories throughout 50 ns simulations run. Glu-170 seemed to produce hydrogen bond in ~51% through water mediated bridges in all the 10416 trajectories. The MD simulations for PKA-Cα- lead1 complex showed similar bonding network with dock complex during all trajectory analysis.
MD simulations clarified the molecular dynamic properties of PKA-Cα active site by the influence of lead1. The analysis of RMSD and RMSF showed the lesser internal motions of side chain residues of receptor PKA-Cα in forming a stable complex with lead1. The optimized potential energy of the PKA-Cα-lead1 system with similar hydrogen bond networking throughout 50 ns simulations run revealed the stability of the complex in a solvated model system. The key residues of allosteric site residues along with the active site residues of PKACα existing between the two domains terminally will be kept opened in forming complex with lead1 (Figure 9). Thus, increased availability of active site and allosteric site residues for the PKA-Cα regulators as a result functioning of PKA-Cα can be elevated. Thus the proposed three leads showed favorable ADME properties and can be used as an agonist to elevate the PKA-Cα expression. Thus elevating the PKA-Cα expression reduces the increased 3R-tau levels results in maintaining a balanced 3R-tau to 4R-tau ratio in the hippocampal neurons ensues the reduced formation of taupathies and neurofibrillary tangles.
For normal functioning of brain the ratio of 3R-tau to 4R-tau should be balanced that is to be 2:1 ratio, but PKA-Cα deregulates the tau splicing alternatively. Thus by elevating the PKA-Cα levels may result in reducing the imbalance in 3R-tau/4R-tau ratio, which ensues the reduced neuronal cell loss with the taupathies. The availability of the multiple co-crystal structures of PKA-Cα was utilized in designing multiple e-pharmacophores based on the structural and interaction energy followed by multiple docking to confess the agonists. Hence in the present study, structure based e-pharmacophore modeling was employed for nine co-crystal structures of PKA-Cα with various filters that included pharmacophore fitness score, XPG Score and ΔG score had conceded potentially binding hits with required ADME/T properties and structural diversity extending to distinguish the novel potential leads as agonists.
E-Pharmacophore modeling has advantages in expeditiously defining the potent hits with structural requirements that are complementary to the site points of the PKA-Cα. The multiple e-pharmacophores demonstrates the different binding modes of the active site residues and allosteric site residues of PKA-Cα and were used in the virtual screening. Multiple e-pharmacophore models provided a frame work for the development of novel leads as an activator by exploiting the key residues from the allosteric site that is outside the binding cleft. The analysis of the lead1- PKA-Cα complexes indicates the key amino acids of allosteric site along with active site of PKA-Cα makes the cleft present between the two terminal domains opened. This result in the elevation of increased availability of active site and allosteric site residues for the regulators of PKA-Cα as a result functioning of PKA-Cα was elevated which results in maintaining a balance between 3R-tau and 4R-tau ratio. Based on RRD, QPLD, IFD and MD simulations results analysis the proposed three leads were enough to elevate the biological activity of PKA-Cα. As PKA-Cα plays a crucial role in maintaining the balance between the 3R-tau and 4R-tau which is necessary for normal brain functioning. PKA-Cα down regulation results in microtubule destabilization and taupathies resulted from increased 3R-tau leading to the loss of hippocampal neurons results in enhancing the memory deficiency in AD. Thus the results emphasized that, the proposed three leads can act as a agonist and also provide a scaffold framework for the elevation of PKA-Cα levels which would be useful for treating PKA-Cα mediated Alzheimer’s disease.
PN is highly thankful to ICMR for sanction of senior research fellowship (BIC/11(04)/2014). Authors are thankful to DBT, Ministry of Science and Technology for providing infrastructure through BIF program (BT/BI/25/001/2006).