Pharmacophore modeling, docking and molecular dynamics simulation for identification of novel human protein kinase C beta (PKCβ) inhibitors.

Protein kinase Cβ (PKCβ) is considered as an attractive molecular target for the treatment of COVID-19-related acute respiratory distress syndrome (ARDS). Several classes of inhibitors have been already identified. In this article, we developed and validated ligand-based PKCβ pharmacophore models based on the chemical structures of the known inhibitors. The most accurate pharmacophore model, which correctly predicted more than 70% active compounds of test set, included three aromatic pharmacophore features without vectors, one hydrogen bond acceptor pharmacophore feature, one hydrophobic pharmacophore feature and 158 excluded volumes. This pharmacophore model was used for virtual screening of compound collection in order to identify novel potent PKCβ inhibitors. Also, molecular docking of compound collection was performed and 28 compounds which were selected simultaneously by two approaches as top-scored were proposed for further biological research. Supplementary Information: The online version contains supplementary material available at 10.1007/s11224-022-02075-y.

Introduction

Protein kinase C (PKC) is a family of serine/threonine protein kinases which transduce essential signaling pathways associated with cell proliferation, differentiation, survival, migration, apoptosis, etc. [1-3]. A total of 15 isozymes of PKC have been reported which are classified based on their structure and mode of activation into 3 subfamilies: classical (PKCα, PKCβI, PKCβII, and PKCγ (activated by diacylglycerol (DAG) and calcium)), novel (PKCδ, PKCε, PKCθ, and PKCη (activated only by DAG)) and atypical (PKCζ and PKCι/λ (activated by protein–protein interactions)) [4]. PKC isozymes are considered as promising molecular targets for cardiovascular [5, 6], neurodegenerative [7, 8], immune [9-11], metabolic [12-14] diseases and different types of cancer [15-17]. Recently, PKCβ has been also considered as a promising molecular target for the treatment of COVID-19–related acute respiratory distress syndrome (ARDS). It was shown that PKCβ is activated by induction of neutrophil receptors (a number of toll-like receptors, CD-14, etc.), and is involved in signaling pathways of neutrophil extracellular traps (NETs) formation. It should be noted that PKCβ activation is specific in induction of NETosis [18-20]. In addition, PKCβ is chronically activated in diabetes due to increased concentrations of DAG [21] and involved in the pathogenesis of diabetic cardiomyopathy [22], diabetic wounds [23], diabetic nephropathy [24] and other diabetic complications [25]. It should be noted that in general, people with diabetes are more likely to have more severe symptoms of COVID-19 [26]. Therefore, an effective inhibition of PKCβ is a promising approach to treat people infected with SARS-CoV-2, especially patients with diabetes.

Nowadays, the inhibitors of PKC have been already published among several chemical classes such as indolocarbazoles [27-30], bisindolylmaleimide derivatives [31-34], balanol analogs [35-38], melittin [39-41], etc. To date, only one selective PKCβ inhibitor, ruboxistaurin which belongs to bisindolylmaleimide derivatives is proposed as investigational drug for diabetic retinopathy. The aim of this study was to perform virtual screening using pharmacophore modeling, molecular docking and molecular dynamics approaches for the search of novel potential PKCβ inhibitors.

Methods

Building of ligand-based PKCβ pharmacophore models

In order to build and validate ligand-based pharmacophore models for human protein kinase C beta inhibitors, the chemical structures of known (PKCβ) inhibitors from the ChEMBL database [42, 43] and literature data [27-41] were collected. The IC50 or Ki values of collected PKCβ inhibitors are available as Supplementary information (Table S1). The structures of the compounds for the training and test sets were minimized (2500 steps) and protonated with OpenBabel 2.4.0 [44] in MMFF94 (mmff94) force field using conjugate gradient algorithm.

At the first step, the primary pharmacophore model of PKCβ inhibitors was constructed based on the intermolecular bonds between bisindolylmaleimide and PKCβ (PDB ID: 2I0E) [45] with Discovery Studio Visualizer 4.0 [46]. PKCβ in this crystal structure has three phosphorylated amino acid residues — Thr500 (in the activation loop of the kinase domain), Thr641 (in the turn motif), and Ser660 (in the hydrophobic region) indicating that the kinase is in an activated form. To calculate pharmacophore feature weights of obtained pharmacophore model, ligand-based pharmacophore modeling using web-server PharmaGist was performed [47]. Bisindolylmaleimide inhibitor, extracted from crystal structure of PKCβ with PDB accession code: 2I0E, which served as a pivot molecule, and 13 the most active inhibitors from different chemical classes with IC50 values less than 5 nM were used for generation of ligand-based pharmacophore models. Other parameters were set by default. Using PharmaGist the superpositions of input ligands with pivot molecule were obtained. The pivot molecule was rigid and the compounds used for superposition were flexible. In this pairwise alignment, all ligands were represented with pharmacophore features. Therefore, the pharmacophore features of the ligands were compared to pharmacophore features of pivot molecule. By default, the distance threshold for the centers of aromatic features was 1.8 Å and for the centers of other types of pharmacophore features was 1.4 Å. The best ligand pairs were selected according to the score value , where — the score of matching for ligand pharmacophore feature with the pivot pharmacophore feature (by default, score value for hydrophobic features is 0.3, for aromatic features — 0.3, for other types of pharmacophore features — 1.0).

Based on the matching of the ligands pharmacophore features to the primary model, the weights of these pharmacophore features were established. The weight of pharmacophore feature means the number of ligands possessing this pharmacophore feature. This parameter indicates the contribution of pharmacophore feature to score during screening.

The models were visualized in Discovery Studio Visualizer 4.0 and saved in .CHM format. Using our in-house program PharmDeveloper [48, 49], the files were converted from .CHM format to .QUERY format. The excluded volumes were added with Discovery Studio Visualizer 4.0. PharmDeveloper program was used for optimization and validation of pharmacophore models, converting of screening database, combining of pharmacophore features, performing of screening and rescoring.

Obtaining of derivative pharmacophore models based on combining of pharmacophore features using PharmDeveloper

If a primary pharmacophore model with n pharmacophore features is given, and a number of possible features in derivative models m is given, then we will get of derivative models with different combinations of pharmacophore features according to formula (1):

The molecule of bisindolylmaleimide, which was used for pharmacophore model construction, was imaginary divided into four substructures — maleimide fragment, aliphatic chain with the terminal amino group and two indoles. The pharmacophore features located within these fragments were united into groups for combining. Each group included a number of pharmacophore features, which were taken for combination. The following parameters for pharmacophore features combination were specified in PharmDeveloper: (a) pharmacophore feature required by default, which will be present in all derived models; (b) feature belonging to a certain group and how many features from this group can be obtained (for example, if there are 3 features in a group and the algorithm can select 2 features from them, we will get 3 different combinations of 2 features from this group ()); (c) features without indicated parameters. The number of required features or required without parameters can be different, including 0. Also, there can be an unlimited number of pharmacophore feature groups (including 0) with a different number of features (at least 1).

First, the algorithm analyzes the primary model. The required pharmacophore features are selected separately and added to all derived models. Each group of pharmacophore features is also analyzed separately and a subset of combinations of pharmacophore features for each group is generated. The number of pharmacophore features obtained from the group of pharmacophore features without parameters is automatically calculated according to formula (2):where is the number of pharmacophore features obtained from the group of pharmacophore features without parameters, — the total number of pharmacophore features in the derived pharmacophore model, — the number of required pharmacophore features by default, — the number of pharmacophore features, obtained from group i.

The next step is the assembly of derived pharmacophore models. These derived models are elements of a multidimensional matrix, each dimension of which is formed by subsets of combinations of pharmacophore features. A special algorithm collects elements of the matrix and adds required by default pharmacophore features to each of them. The number of obtained derived models is calculated according to formula (3):where is the number of derived pharmacophore models, — number of combinations, obtained from features without parameters, — number of combinations, obtained from group i.

Optimization of ligand-based PKCβ pharmacophore models

Optimization is an iterative procedure in which each step of the iteration is a validation screening of the test set toward optimized model. We have performed three iterations of pharmacophore model score and radius optimization with the steps of 0.1 and 0.2 for each parameter, respectively. For pharmacophore model optimization, we used chemical structures of reported protein kinase Cβ (PKCβ) inhibitors obtained from ChEMBL database and literature data. Totally there were 731 compounds. The training 13 inhibitors, which were used during pharmacophore model construction, were removed from this set. The remaining compounds were divided into active and inactive. The cutoff of activity was ≤ 50 nM. There were 303 active compounds and 415 non-active compounds.

During the validation of pharmacophore models, the ligands were evaluated by the PDscore value, which demonstrates the correspondence of the ligand to the pharmacophore features of the model. The more precisely the ligand corresponds to all pharmacophore features — the score is higher. Besides geometric matching, also the weights of pharmacophore features are taken into account. Matching with a feature of greater weight gives a greater increase in score than matching with a feature of smaller weight. Also, an important component of this score is the QSAR assessment based on active/inactive ligands from the test set (for this purpose, the program creates the table of special reference coefficients to which the parameters of the investigated ligands are compared). This process is fully automated. PDscore was calculated by the formula presented in the Supplementary information. For QSAR assessment molecular descriptors such as number of atoms, number of bonds of different types, molecular weight, number of donors and acceptors of hydrogen bonds, logP, number of halogen atoms, topological polar surface area (TPSA) and molecular refraction were calculated for each molecule.

Molecular docking

The semi-flexible molecular docking was performed with Autodock 4.2.6 [50]. The crystal structure of human protein kinase Cβ (PKCβ) was obtained from the Brookhaven Protein Data Bank (PDB ID: 2I0E) [45]. The ligand – bisindolylmaleimide inhibitor, water molecules and ions were removed from the PDB file using Discovery Studio Visualizer 4.0 [47]. The receptor was prepared with MGL Tools 1.5.6 and AutoGrid within Autodock 4.2.6 software suite. Ligands were prepared by OpenBabel [44], Vega ZZ (command line) [51] and MGL Tools 1.5.6. The incoming formats of receptor and ligand molecules were converted into PDBQT formats with Vega ZZ using AUTODOCK force field. For molecular docking, parameter files were created for each ligand in the DPF (Dock Parameter File) format, where all necessary parameters for calculations were specified. These files were prepared by scripts included in the MGL Tools software package. The parameters for virtual screening were the following: translation step — 2 Å, quaternion step — 50°, torsion step — 50°, torsional degrees of freedom and coefficient — 2/0.274, cluster tolerance — 2 Å, external grid energy — 1000, max initial energy — 0, max number of retries — 10,000, number of individuals in population — 300, maximum number of energy evaluations — 850,000, maximum number of generations — 27,000, number of top individuals, which survived to the next generation — 1, rate of gene mutation — 0.02, rate of crossover — 0.8, mode of crossover — arithmetic. Alpha parameter of Cauchy distribution was 0, Beta parameter Cauchy distribution — 1. The number of iterations of the Lamarckian genetic algorithm was 10 for each ligand. The ranking of docking results was performed by score values calculated by the AutoDock scoring function. The best-scored complexes were used for visual inspection.

Visual analysis of the docking results was done using Discovery Studio Visualizer 4.0 [47]. The complexes were evaluated by the ability of ligands to form hydrogen bonds with the key amino acid residues in the ATP-binding site of PKCβ.

Molecular dynamics simulation

Molecular dynamic (MD) simulations were performed using GROMACS 4.5 software package [52-54]. The starting coordinates were taken from the crystal structure of human protein kinase PKCβ with bisindolylmaleimide inhibitor (PDB ID: 2I0E). Topology file for ligand was obtained with Dundee PRODRG server [55]. We used GRO file of the ligand containing polar/aromatic H’s. Topology file for PKCβ was generated from PDB-file using pdb2gmx command. MD simulations of receptor-ligand complex with explicit inclusion of water were carried out using the GROMOS96 force field [56]. Energy minimization of receptor-ligand complex was done with steepest descent algorithm for 1000 relaxation steps. Then, position restrain dynamics of minimized structure was performed for 20 ps. After that, MD simulation of the system was carried out for 20 ns at constant volume and temperature (NVT). The integration of the equations of motion was done with the standard leap-frog algorithm [57]. The Particle Mesh Ewald (PME) algorithm [58, 59] was used for electrostatic processing.

Results and discussion

The primary pharmacophore model of PKCβ inhibitors was constructed based on the intermolecular interactions of bisindolylmaleimide inhibitor with amino acid residues in the ATP-binding site of human PKCβ (crystal structure with PDB accession code: 2I0E) using Discovery Studio Visualizer 4.0. The interactions of the ligand with amino acid residues of PKCβ active site are shown in Fig. 1. All aromatic rings of the two indole and one maleimide fragments of the bisindolylmaleimide inhibitor form a number of hydrophobic interactions with amino acid residues Ala369, Val356, Leu348, Val423, and Ala483 in the ATP-binding site of protein kinase PKCβ. These interactions were represented by five aromatic pharmacophore features without vectors with the centers, located in geometric centers of inhibitor aromatic rings with standard radius of 1.1 Å. Methyl at the C2-position of indole ring forms hydrophobic interactions with Met473 and Ala483, therefore we built hydrophobic pharmacophore feature on this substituent with radius of 1 Å. Two hydrogen bond acceptor pharamacophore features were generated on two keto groups of maleimide which form hydrogen bonds with Val423 and Thr404. Two hydrogen bond donor pharamacophore features were built on nitrogen atom of maleimide which forms hydrogen bond with Glu421 and on the terminal amine group of inhibitor which is involved in hydrogen bond formation with Asp470. The radii of hydrogen bond acceptor and donor pharamacophore features were 0.91 Å (Fig. 1). Also, excluded volumes were added on all atoms of amino acid residues of PKCβ ATP-binding site in the radius of 5 Å around the ligand. The excluded volumes make the model more accurate since they restrict selection of compounds which will potentially overlap with the amino acid residues in the active site. Hence, the selected compounds will correspond to the volume and shape of ATP-binding pocket. Actually, the shell of the excluded volumes mimics the active site of the receptor. The radius of the excluded volumes was set to 1.2 Å, which corresponds to the van der Waals radius of the hydrogen atom.

Fig. 1

The complex of bisindolylmaleimide inhibitor with amino acid residues in the ATP-binding site of PKCβ. The hydrogen bonds are represented by green dashed lines and hydrophobic interactions are shown by magenta dashed lines

As a result, the primary PKCβ pharmacophore model consisted of ten pharmacophore features including five aromatic features without vectors, one hydrophobic feature, two acceptors and two donors of hydrogen bonds and 158 excluded volumes (Fig. 2). In order to make presentation of pharmacophore model clear, excluded volumes are not visualized since they cover the pharmacophore features.

Fig. 2

The primary pharmacophore model of PKCβ inhibitors. Aromatic pharmacophore features (Ar1-5) without vectors are indicated with blue color, hydrophobic pharmacophore feature (Hyd) is presented by cyan color, hydrogen bond acceptor pharmacophore features (Acc1-2) are shown by green color, hydrogen bond donor pharmacophore features (Don1-2) are labeled with magenta color

The primary model had default radiuses of pharmacophore features and did not have weights for pharmacophore features. In order to calculate this parameter, we performed ligand-based pharmacophore modeling using web-server PharmaGist. For this, bisindolylmaleimide inhibitor extracted from crystal structure with PDB ID: 2I0E (served as a pivot molecule) and 13 of the most active compounds of this protein kinase with IC50 value less than 5 nM, which were obtained from ChEMBL database, were uploaded in .MOL2 format and submitted in PharmaGist. The pairwise superpositions of the ligands with bisindolylmaleimide were performed and the weights of pharmacophore features were determined. The weight of pharmacophore feature corresponds to the number of ligands possessing this pharmacophore feature. This parameter means the contribution of pharmacophore feature to score value during screening. The pharmacophore features weights are presented in Table 1.

Table 1

Pharmacophore feature weights of the primary PKCβ inhibitors pharmacophore model

Pharmacophore feature	Acc1	Acc2	Don1	Don2	Ar1	Ar2	Ar3	Ar4	Ar5	Hyd
Weight	7	7	9	5	9	7	4	9	6	3

Since the primary model was very complex — it consisted of ten pharmacophore features and 158 excluded volumes, we obtained simpler derived pharmacophore models using algorithm of pharmacophore features combining in PharmDeveloper program. The total number of pharmacophore features in derived models was set as five or six. During the combining procedure, some pharmacophore features were set as required by default, several pharmacophore features formed combination group, from which some features were selected and others were removed. There were 9 different directions of combining (Table 2). Totally, 208 derived different pharmacophore models were generated.

Table 2

The combining directions of pharmacophore features for obtaining of derived pharmacophore models

Combining direction	1	2	3	4	5	6	7	8	9
Number of pharmacophore features in derived models	5	5	6	6	5	5	5	5	5
Number of derived models	24	24	24	24	24	24	16	24	24
Pharmacophore features and its parameters for combination
Acc1	g1-1	g1-1	g1-2	g1-2	g1-2	g1-2	g1-1	g1-2	g1-2
Don1	g1-1	g1-1	g1-2	g1-2	g1-2	g1-2	n	g1-2	g1-2
Acc2	g4-1	g1-1	g1-2	g1-2	g1-2	g1-2	g1-1	g1-2	g1-2
Ar1	n	n	n	n	-	-	n	n	n (Hyd)
Ar2	g2-1	g2-1	g2-1	g2-1	g2-1	g2-1	g2-1	g2-1 (Hyd)	g2-1 (Hyd)
Ar3	g2-1	g2-1	g2-1	g2-1	g2-1	g2-1	g2-1	g2-1 (Hyd)	g2-1 (Hyd)
Ar4	g3-1	g3-1	g3-1	g3-1	g3-1	g3-1	g3-1	g3-1 (Hyd)	g3-1 (Hyd)
Ar5	g3-1	g3-1	g3-1	g3-1	g3-1	g3-1	g3-1	g3-1 (Hyd)	g3-1 (Hyd)
Hyd	g4-1	g4-1	g4-1	g1-2	g4-1	g1-2	g1-1	g1-2	g1-2
Don2	g4-1	g4-1	g4-1	n	g4-1	n	g1-1	-	-

gX-Y, combination group with name X, from which algorithm will select number Y of pharmacophore features in each iteration of combination; n, is feature required by default, which will be present in all derived models; (Hyd), aromatic feature is changed to hydrophobic feature

These obtained derived models were validated on a test set, which included all known PKCβ inhibitors collected by us from ChEMBL database and the articles [27-41], except the most active 13 compounds which were used as a training set during primary pharmacophore model construction. These ligands were divided into active and inactive using the activity cutoff of ≤ 50 nM. Therefore, there were 303 active compounds and 415 non-active. All 208 derived models were used for pharmacophore screenings of the test set. It was found that only 19 models correctly identified at least one or more active compounds. Then, these 19 models were optimized and validated using the same test set. The optimization is an iterative procedure, in which each iteration step is a validation screening of the test set against the optimized model. Three iterations of radius and score optimization were performed with step = 0.1 and 0.2, for each parameter, respectively. Therefore, as a result we obtained 1026 derivative models with different parameters. Among them, we selected 62 models which correctly identified more than 60% active compounds in the test set. The parameters of the models are available as Supplementary information (Table S2). The most accurate PKCβ pharmacophore model, which correctly predicted more than 70% active compounds of test set (in other words, the model had the percentage of active compounds among all selected ligands (including non-active) more than 70%) is presented in Fig. 3. Mapping of the 13 most active PKCβ inhibitors that were used as the training set, 10 active and 10 non-active PKCβ inhibitors from the test set to final pharmacophore model are available as Supplementary information (Figs. S1, S2, S3).

Fig. 3

The validated ligand-based pharmacophore models of PKCβ inhibitors (PKCb_2-5combf-22_624_m (a), PKCb_2-5combf-23_372 (b)). Aromatic pharmacophore features without vectors are indicated with blue color, hydrophobic pharmacophore feature is presented by cyan color, hydrogen bond acceptor pharmacophore features are shown by green color, hydrogen bond donor pharmacophore features are labeled with magenta color and excluded volumes are not shown since they cover all the model

In order to establish whether pharmacophore models are stable during protein kinase C beta conformation changes, we have performed molecular dynamics simulation of the complex of human PKCβ with bisindolylmaleimide inhibitor during 20 ns. The RMSD plot of the ligand, sum of Coulomb and Lennard–Jones interaction energies for inhibitor-PKC complex and hydrogen hydrogen-bond existence map are shown in the Fig. 4. The presented hydrogen bond network (Fig. 4c) correlates well with the energy of the complex (Fig. 4b). The decrease of the energy of the complex, which means further stabilization, can be associated with the disruption of hydrogen bond of the ligand with amino acid residue Val423 (Fig. 4c, E) and by formation of hydrogen bonds with Tyr422 (Fig. 4c, D) and Arg624 (Fig. 4c, F-L). Then, the breaking of these hydrogen bonds at 4.5 ns leads to increase of the energy but further formation with Tyr422 at 6 ns again causes decrease of the energy. As it can be seen from the RMSD plot, some disposition of compound in the ATP-binding site is observed at 14 ns of MD simulation (Fig. 4a). At this time, the stabilization of PKC-inhibitor complex obviously is supported by the formation of hydrogen bonds with Lys350 (Fig. 4c, A, R), Lys371 (Fig. 4c, B), Asp484 (Fig. 4c, O, P) and Leu348 (Fig. 4c, Q). It should be noted, that the inhibitor forms hydrogen bonds with Asp484 (Fig. 4c, Y, Z) during almost all the time of MD simulation, except the range from 16 to 19 ns, where these hydrogen bonds are substituted by hydrogen bonds with Asp470 (Fig. 4c, T, U). Therefore, hydrogen bonds of the ligand with hydroxyl group of Tyr422 and carboxyl group of Asp484 are very important for PKC-inhibitor complex stabilization.

Fig. 4

Molecular dynamic (MD) results of bisindolylmaleimide inhibitor in complex with PKC during 20-ns MD simulation: RMSD of inhibitor (a), sum of Coulomb and Lennard–Jones interaction energies for inhibitor-PKC complex (b), hydrogen-bond existence map of inhibitor-PKC complex: A – Lys350_N – O14, B – Lys371_NZ – O20, C – Tyr422_OH – N13, D – Tyr422_OH – O14, E – Val423_N – O14, F – Arg624_NE – O20, G – Arg624_NE – N13, H – Arg624_NH1 – O20, I – Arg624_NH1 – N13, J – Arg624_NH2 – O20, K – Arg624_NH2 – N13, L – Arg624_NH2 – O14, M – Asp470_OD1 – N25, N – Asp470_OD2 – N25, O – Asp484_OD1 – N25, P – Asp484_OD2 – N25, Q – Leu348_O – N13, R – Lys350_O – N13, S – Asp470_OD1 – N30, T – Asp470_OD2 – N30, U – Asp470_O – N30, W – Asn471_OD1 – N30, V – Asn471_ND2 – N30, X – Ala483_O – N30, Y – Asp484_OD1 – N30, Z – Asp484_OD2 – N30 (c)

The MD frame of PKC-inhibitor complex with the lowest energy value at 4000 ps of MD simulation was used to check how it corresponds to the best pharmacophore models. The superposition of the starting ligand conformation, extracted from crystal structure (PDB ID: 2I0E), which was used for pharmacophore models construction, and ligand conformation obtained from minimized complex, demonstrates that the ligand geometry changed minimally (RMSD = 0.33 Å) after MD simulation (Fig. 5a), while the interactions with amino acid residues in the active site of PKCβ changed significantly. As it can be seen from Fig. 5b, these changes did not impact the model PKCb_2-5combf-23_372, since the part of the ligand, where pharmacophore features were added, is very stable. Instead, the model PKCb_2-5combf-22_624_m was changed due to displacement of the flexible terminal amino group of the inhibitor. The hydrogen bond donors which were changed are shown by the gray pharmacophore features (Fig. 5c). Therefore, we have selected the model PKCb_2-5combf-23_372 for pharmacophore screening since it is very conservative during conformation changes of protein kinase and has a good hit rate in experiments on the test set.

Fig. 5

The superposition of the starting ligand conformation, extracted from crystal structure (PDB ID: 2I0E), which was used for pharmacophore models construction (carbon atoms are labeled with red color), and ligand conformation obtained from minimized complex (carbon atoms are labeled with green color) (a), mapping of the ligands to the model PKCb_2-5combf-23_372 (b), mapping of the ligands to the model PKCb_2-5combf-22_624_m (c). Aromatic pharmacophore features without vectors are indicated with blue color, hydrophobic pharmacophore feature is presented by cyan color, hydrogen bond acceptor pharmacophore features are shown by green color, hydrogen bond donor pharmacophore features are labeled with magenta color and the hydrogen bond pharmacophore features which were changed are shown by the gray color

The obtained pharmacophore model which included three aromatic pharmacophore features without vectors, one hydrogen bond acceptor pharmacophore feature, one hydrophobic pharmacophore feature and 158 excluded volumes, was used for virtual screening of the commercially available OTAVA compound collection [60] containing about 150,000 compounds in order to select the compounds for study of inhibitory activity toward PKCβ. According to the values of RMSD, PDscore and visual analysis of compounds for matching with pharmacophore model, we have selected 4120 compounds.

Also, we have performed molecular docking of OTAVA compound library into ATP-binding pocket of PKCβ. According to score values and ligand-receptor intermolecular interactions such as formation of hydrogen bonds with amino acid residues in the hinge region and hydrophobic interactions with Leu348, Ala369, Ala483, we have selected 651 compounds. These compounds were compared with the top-ranked compounds selected according to pharmacophore screening. It was found that both sets of compounds include 28 identical molecules which can be proposed for biochemical screening. The chemical structures of compounds are presented in Table 3.

Table 3

The structures of top-scored compounds selected simultaneously by pharmacophore screening and molecular docking

Conclusion

The ligand-based pharmacophore model of human PKCβ was built and used for virtual screening of OTAVA compound collection. Also, the molecular docking was performed. The compounds which were selected simultaneously by two approaches as top-scored were proposed for further biological research.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 934 KB)