In Silico Insights of Tunable Small Molecules With Regulatory Activity of Acid Sphingomyelinase- Protein Kinase Cδ Interaction


 Acid sphingomyelinase demonstrates a housekeeping role and plays a key role in maintaining membrane turnover through sphingolipid metabolism. While deficiency of ASM housekeeping function leading to a lysosomal storage disorder, activation and translocation of ASM leads to AD pathology. Phosphorylation of ASM by PKCδ is critical for translocation of ASM to the plasma membrane. However, several strategies have been developed for the treatment of patients with AD, limitations of treatments and occurrence of AD by alternative signaling pathways become hurdle to the mortality of patients. In the present study, we report the interface of ASM- PKCδ interactions by performing molecular modeling, docking and dynamics simulation analysis. By considering the interacting site of ASM with PKCδ as an allosteric site, we screened a large number of small molecules by virtual screening and identified four lead molecules that show high binding affinity and interacting mode with residues of ASM that are critical to making interaction with PKCδ. The proven role of ASM in apoptosis, neurodegenerative diseases, and the existing functionality of leads in these diseases, together with strengthens the lead compounds, ZINC85551993, ZINC71928291, ZINC71773625 and ZINC85432419 as allosteric inhibitors for ASM. We believe that the lead molecules could be potential allosteric modulators for the translocation of ASM.


Introduction
A balanced Sphingolipid metabolism is essential to accomplish numerous cellular processes for the maintenance of normal physiological conditions. However, dysregulation of sphingolipid cues leads to numerous diseases including neurodegenerative complications. Alzheimer's disease (AD) is the most widespread dementia, accounting for up to 80% of dementia cases and it is more prevalence in aging population [1] [2]. AD is characterized by the accumulation of misfolded proteins, including amyloid β, thereby neuro brillary tangles [3] [4]. While post-mortem evaluation is the de nitive diagnosis of AD, recent research advances the diagnosis of AD in living patients by the detection of positron emission tomography (PET) and cerebrospinal uid (CSF) biomarkers [5]. Additionally, more recently, several approaches have been established for diagnosis of AD, including non-invasive imaging to detect amyloidβ, evaluation of hyperphosphorylated tae peptide (p-tae), assay for circulating proteins of serum, and serum microRNA pro ling [6] [7] [8] [9]. Along with establishment of diagnosis, researchers have been developed therapeutic approaches for the treatment of patients with AD. Cholinesterase inhibitors, including galantamine, donepezil and rivastigmine are currently recommended for the treatment of patients with mild-to-severe AD [10]. Memantine, an antagonist of N-methyl-D-aspartate receptor and agonist of dopamine is approved for the treatment of patients with AD [11]. In addition, huperzine A, omega fatty acid supplementation are also received attention in the implementation of AD treatment strategies [12] [13] [14].
Various lipid classes are much abundant in the cell membrane composition, and along with protein/receptors of cell membrane, these lipids demonstrate numerous cell signaling events to coordinate cellular organelles and cell environment [15]. Perturbation of membrane lipids leads to adverse effects and is the major underlie cause of AD pathology. Ceramide, a class of lipid has received greater attention among the lipid classes as it contributes to AD pathology by affecting amyloid-β generation and phosphorylation [16] [17]. Emerging evidence has revealed that the lowering ceramide levels could reduce accumulation of amyloid-β, and indicated the importance of ceramide in therapeutic effectiveness [18] [19]. As a central signaling molecule in sphingolipid metabolism, ceramide can be produced by degradation of sphingomyelin, or synthesized de novo from sphingosine and fatty acyl coA [20] [21] [22]. Existing evidence has reported that acid sphingomyelinase (ASM) that involved in sphingomyelin hydrolysis, and acid ceramidase that responsible for ceramide synthesis, are expressed at elevated levels in AD [16] [23]. ASM, localized to lysosomes with an optimal pH, is play key role in ceramide formation by degradation sphingomyelin [24]. Moreover, an increased ceramide formation in lipid raft microdomins by the induction of activation and translocation of ASM has been identi ed [25]. Protein kinases are well-known to their indispensable role in activating substrate proteins. Accumulating evidence has revealed that protein kinase C delta (PKCδ) activates and translocate to plasma membrane by phosphorylation [26]. While the overexpression of PKCδ associated with neuronal damage, its inhibition results in reduction of cellular injury has been reported [27] [28] [29].
Although, several approaches has been developed in recent years to treat patients with AD, ine cacy of disease modifying and limitations of treatment, are still hurdle in mortality of patients with AD [30].
Moreover, identifying the potential signaling pathways in AD pathology is being increased for the implementation of AD clinical settings. Based on the facts, we believe that ASM-PKCδ signaling is a major contribution to AD pathology. In this study, we explore ASM-PKCδ binding interface by molecular docking and molecular dynamics simulation. We further screened a large number of small molecules against PKCδ binding site of ASM to identify lead molecules as the allosteric inhibitors.

Materials And Methods
Full-length Structure prediction of PKCδ and ASM In order to understand the molecular interaction between PKCδ and ASM, the full length model of PKCδ was predicted by using Robetta server [31] [32] that demonstrates homology modeling and ab initio methods. The sequence of PKCδ (UniProt ID: F1MZX0, 676aa) was retrieved from the UniProt database (http://www.uniprot.org) and subsequently performed protein BLAST against the Protein Data Bank (PDB, http://www.rscb.org/pdb/) to identify homologue templates. Due to lack of homologous template, Robetta server was used to model the protein structure. Robetta server aligns PKCδ sequence with extended-synaptotagmin (PDB ID 4P42), protein kinase C γ (PDB ID 2E73), human chymerin1 (PDB ID 3CXL), lambda cro (PDB ID 3A63) and protein kinase alpha (PDB ID 3IW4). Based on the essentiality of zinc metal ion for its activity, the obtained model was remodeled using Modeler9v8 program [33] to x four zinc metal ions. As a nal step of structure prediction, the least Modeler Objective Function (MOF) and Discrete Optimized Protein Energy (DOPE) [34] score were considered to select the best model selection. Further, the model was validated by the Structural Analysis and Veri cation Server (SAVES) server that contains PROCHECK and ERRAT [35].
The full length ASM structure prediction was clearly described in our previous study [36]. However, recently the crystal structure of human acid sphingomyelinase (PDB ID 5JG8) with unsolved small portions of N-terminal and C-terminal regions was deposited in the protein data bank [37]. Thus, the crystal structure was taken into consideration as a template to remodel our previous ASM model using Modeler9v8 program. The aforementioned procedure was followed for the selection of the best model, and subsequent structural veri cation.

Structural preparation of PKCδ
To attain active state of PKCδ, docking was performed with ATP molecule at the region of ATP binding site; ATP (CID: 5957) was downloaded from PubChem database (http://pubchem.ncbi.nlm.nih.gov/search/search.cgi) and followed ligand preparation using LigPrep module of Maestro (Schrödinger 2014-2) software by assigning OPLS-2005 force eld. At this point, the most occurring possible ionization states at pH7.0 ± 2.0 were set up with their speci ed chiralities. The PKCδ protein was prepared to adapt the bond order and setting the hydrogens using the protein preparation wizard, which is embedded in Maestro. Further, hydrogen bonds were assigned based on the sample water orientation, and optimized by exhaustive sampling methods. Moreover, the protein was subjected to energy minimization by applying OPLS 2005 force eld. Residues (355-363) of PKCδ, those proven as ATP interacting residues were selected to generate receptor grids. The grid of about 10 Å radius was created to dock the molecule by following the standard precision (SP) protocol of Maestro.

Molecular dynamics simulation of PKCδ and ASM
To check the structural stability, MD simulation of PKCδ and ASM was performed using GROMACS 5.1.0 molecular dynamics package [38] [39]. The topology information of ATP molecule was generated based on PRODRG Server and the topology les of PKCδ and ASM were generated using pdb2gmx command. GROMOS96 54a7 force eld was used to assign the proper bond order and geometrical parameters of the proteins. Moreover, the systems with cubic box (distance of 1.0 nm from the center of the protein to edge of the box) were solvated using SPC216 water model, and periodic boundary conditions were applied to maintain the directions of systems. Based on the net charge of systems, water molecules were replaced with Cl-ad Na + counter ions. Next, energy minimization of proteins was performed using steepest descent algorithm with a tolerance of 1000kJmol − 1 nm − 1 . Moreover, interactions including electrostatic and Van der Waals were applied by using Fast particle-Mesh Ewald electrostatics (PME) [40] with 1nm cut off. LINCS and SETTLE algorithm were used to constrain the bond angles and geometry of the water molecules, respectively. While, Parrinello-Rahman [41] method was used regulate the pressure, the modi ed weak coupling berendsen thermostat and V-rescale algorithm were used to regulate the temperature of the systems. The distance restraints method (the distance ranging between 2.3 and 3.3 Å) was used to retain the coordination bonds formed between Zn2 + metals and metal binding residues. The residues including D202, H204, D274, N314, H421, H453 and H455 of ASM, and C189, C208, C192, H159, C172, C200, C175, C280, C264, C261, H231, C244, C247 and C272 of PKCδ were involved in the formation of coordination bonds. Both NVT (constant number of particles, volume and temperature) and NPT (constant number of particles, pressure and temperature) were accomplished for 100 ps, and checked for their equilibration status. After equilibration, MD run for about 50ns of production was performed with a time frame of 2 fs. The dynamics and structural transitions of PKCδ and ASM were analyzed using the tool of GROMACS package.

The PCA and FEL analysis
The high amplitude concerted motions of protein in simulated trajectories through eigenvectors of the covariance matrix were determined by using PCA and dPCA analysis [42]. Accordingly, the PCA of PKCδ and ASM proteins were calculated to demonstrate the atomic motions of the proteins during the course of MD simulation. The cosine content (cj) was calculated for each principal component (pi) of covariance matrix to signify the free energy landscape [43] [44]. These sensitive measures of trajectories enable mapping free energy landscape of chosen principal components (PCs), which describes further conformational changes of proteins. In general, the values of cosine content ranging between 0 (no cosine) and 1 (perfect cosine) in the total time of simulation period (T): While, the structural dynamics of protein can be understand by the rst eigenvector cosine contribution, in most of the cases, the value of the rst eigenvectors cosine shows closer to 1. As the high cosine content represents the large-scale motions, the rst cosine content, generally, is not suitable for analyzing FEL.
However, recent studies reported that the cosine contents possess the value ranging below 0.2 or sometime up to 0.5 can result accurate conformational mapping [44]. Therefore, here, the rst 20 PCs of each protein were analyzed based on their cosine values. FEL was prepared using cosine contents, which showed lesser than 0.2 of the rst two projection eigenvectors and were represented as PC1 and PC2 respectively. These analyses provided the energy minima of each protein, which further led to point out the best free energy conformation of PKCδ and ASM.

Protein-Protein interaction studies
To understand the structural interacting pattern of ASM, protein-protein docking was performed between PKCδ with ASM using the High Ambiguity Driven Bio-molecular Docking (HADDOCK) [45] [46]. The most stable conformations of both ASM and PKCδ with lowest free energy were prepared by assigning the hydrogen atoms followed by minimization. PKCδ interacting region of ASM from literature was used as active site residues for docking purpose, while the neighboring amino acids of the active site were taken as passive residues. The best cluster was selected based on the HADDOCK score and RMSD cut-off of 7.5 Å, and further analysis was performed by using PISA (Protein Interfaces, Surfaces and Assemblies) web interface. A binding free energy of docked complex was calculated using PDBePISA interface [47] [48]. Moreover, molecular dynamics simulation of ASM-PKCδ complex was performed to understand the structural transitions of ASM upon interaction with PKCδ.
Free energy calculations Molecular mechanics Poisson Boltzmann or Generalized Born solvent accessible surface area [MMPB (GB) SA] [49] [50] is a computational tool used for the analysis of bio molecular interactions. The principles behind this method are well authorized and used effectively in earlier studies [49]. In recent years, several drug discovery studies have used this tool as a scoring function. Here, the free energy of protein-protein complex was obtained by using MMPBSA algorithm through Gromacs 5.1.0. The equilibration region (50-80 ns) trajectory was used for this application. All water molecules and ions were removed from the trajectory. The free energy calculations were determined by the following equation:

Molecular interactions of ASM allosteric site
To identify allosterically binding lead compounds, lowest energy conformation of ASM was selected through MDS and PCA analysis. Accordingly, the protein was prepared in multistep process using protein preparation wizard of Schrodinger software. After adding hydrogen and bond order assignment, the protein was subjected to energy minimization by OPLS2005 force eld. The grid was generated by using allosteric site information from ASM-PKCδ docking and experimental evidence. A total of 958771 library of compounds from zinc data base were retrieved and prepared via ligprep module in Schrodinger software. The prepared ligand compounds were docked into ASM allosteric site through High Throughput Virtual Screening (HTVS) using glide maestro module. The successful 10% compounds of HTVS were taken up to SP docking for further screening. Finally, to eliminate false positive results, the best 10% compounds from SP docking were allowed for XP docking. The results obtained in virtual screening were analyzed. The gures were prepared using glide maestro module of Schrodinger and chimera software.

Results And Discussion
Structural analysis of PKCδ The modeled full length PKCδ protein indicates that the structure is composed of regulatory (1-329) and catalytic (330-676) subunits. The regulatory subunit consists of three small domains include a C2 domain (1-146) which is believed to be non-calcium binding domain and two Zinc nger motifs (158-208 and 230-280). The C2 domain comprises of 8 β strands that are involved in the formation of β sheet, and one helix. The catalytic subunit is composed of a protein kinase domain (349-603) and an AGCkinase C-terminal domain (604-676) (Fig. 1). The protein kinase domain is composed of N-lobe and Clobe and these two lobes are joined by a hinge linker region. This type of arrangement of protein kinase domain of bovine PKCδ shows homologue nature with other protein kinases. Six-stranded β-sheet and three α-helices collectively forms the N-lobe whereas the C-lobe consists mostly of helical structural elements (8 helices) and two β-strands. The ATP molecule binds in the ATP binding region (355-363) through hydrogen bond formation. Consequently, the β phosphate oxygen atom O8 forms hydrogen bonds with K378 and A490 with contact distance of 2.92 Å, 2.44 Å whereas the O5 oxygen atom forms hydrogen bonds with K378 with a contact distance of 3.04 Å respectively. The O2 and O3 atoms of the ribose moiety forms two non-classical hydrogen bonds with the backbone oxygen atom of N478 with a distance of 2.96 Å and 2.64 Å respectively. The NH2 atom of purine ring forms classical hydrogen bond with D634 oxygen atom of PKCδ with a distance of 3.03 Å. In addition, F633 (CZ) forms T-shaped π-π interaction with Cg of purine and pyrimidine moiety with contact distance of 4.1 Å and 3.8 Å respectively. The residues L401, F360, L480, V363, A376, A490, F637 and M427 forms hydrophobic interactions whereas T411, N478, K475, E397 and D477 forms polar and charged interactions with ATP (Fig. 2). All these interaction analysis shows that ATP resides in perfect binding mode at the ATP binding cavity as like in other PKC kinase family proteins with the Glide energy of -56.975 kcal/mol.

Structural analysis of ASM
The full length ASM model details were described in our previous publication [36]. However, after remodeling based on recently available human ASM crystal structure, the ASM structure has adapted some additional features. The full length ASM is comprised of three different domains, which includes Nterminal domain (NTD), phosphoesterase domain (PED) and C-terminal domain (CTD). The helical domains (contains only helices) NTD and CTD are positioned at the angle and distance of 45.6˚ and 23 Å with respect to the domain axis. The PED domain is placed with angle and distance of 83.2˚ and 7.8 Å from NTD and 50.0˚ and 16 Å from CTD (Fig. 3). The N-terminal domain comprises a total of 6 helices. The helices α1 (C29-A42) and α2 (K66-L74) were placed parallel to each other. The helices α3 (T84-L99), α4 (Q102-L119), α5 (P124-T143) and α6 (a152-L156) contributes to the formation of saposin B type domain and plays a role in lipid binding. Unlike, previous ASM model, saposin B type domain of remodeled ASM adopted open con rmation and this domain has moved towards the catalytic cavity. This clearly indicates that after lipid binding, saposin B type domain handover the lipid molecule to the catalytic site for its hydrolysis. The PED contains totally 9 helices (α 7-α 15) and 9 beta strands (β 1-β 9). Among 9 β strands, β4, β5, β6 and β7 have contributed to the catalytic site and oriented in such a way that the residues H421, H453 and H455 have coordinate interactions with zinc metal ions along with H203, D202, D274 and N314 residues. Particularly, the long loop region (D206-L248) has adopted open conformation by giving a way for the substrate entry into catalytic cavity. This con rmation is similar to ASM structure at pH 5.0 of our previous study.

Molecular dynamics simulation
To understand the structural stability of PKCδ and ASM, protein structures were subjected to molecular dynamics simulation studies at 50 ns for each protein. The backbone RMSD of PKCδ has shown higher uctuations (0.5-1.75nm) upto 25 ns and then it was well equilibrated till 50 ns MD run. The backbone RMSD pro le clearly indicates that PKCδ protein attains stable conformation from 25-50ns production MD run with little RMSD uctuations (Fig. 4A). Further, the RMS uctuation data of PKCδ indicates that more number of residues (125-150) of the C2 domain have higher uctuations upto 1.8 nm whereas in other domain (zinc nger domains and protein kinase domain) only few residues have shown higher uctuations. (Fig. 4B). The same analysis was carried out for ASM to know the protein stability by means of RMSD uctuations. The RMSD pro le of ASM shows that the protein took rst 15 ns time period to reach equilibration state and subsequently equilibrated till 50 ns MD run with very less RMSD uctuations. This clearly depicts that ASM is energetically stable for about 35 ns time scale of MD run (Fig. 4A). The RMS uctuation pro le shows that maximum number of residues of the signal peptide has higher uctuations up to 1.25 nm. Residues of saposin B type domain have uctuation up to 0.5 nm, and the residues of other domains like PED and CTD have very less uctuations compared to NTD (Fig. 4C).
This analysis clearly depicts the overall stability of individual domains of PKCδ and ASM.

Structural transitions of PKCδ
In order to understand the structural transitions attained in the MD run production, representative structures were extracted and explored based on the Principle Component Analysis and Free Energy Landscape graph. Accordingly, the rst two principle components that show less than 0.2 cosine content were used to construct the FEL graphs. The contour map of the FEL analysis showed three different minimum energy clusters (Fig. 5). Representative structure of cluster 1 indicates that the C2 domain (C2D) is placed at about 27.1 Å distance from the Zinc nger domain (ZFD) whereas representative structures of cluster 2 and 3 moved closer with distance of 21.0 and 23.5 Å respectively, based on the domain-domain distance analysis ( Table 1). The ZFD and Catalytic domain (CD) distance in cluster 1, 2 and 3 were about 27.4, 28.3 and 27.4 Å respectively, which clearly shows the compact nature of ZFD with CD throughout the time period of production MD run. In cluster1 and 2 the catalytic domain is placed at about 39.9 and 36.8 Å distance respectively with C2D whereas in cluster3 the C2D moved away from the catalytic domain with a distance and angle of 43.7 Å and 18.9˚ respectively. The AGC kinase domain loses its secondary structural elements in the cluster2 representative structure. Thus, Based on the domain distance analysis of clusters, cluster 1 was found to be best for further docking studies.

Structural transitions of ASM
In the case of ASM structure, FEL analysis showed two free energy minimum clusters based on the rst two principle components. The representative structure of each free energy minimum cluster was extracted and displayed in Fig. 6. NTD maintained an average distance of 24.1Å from PED in both FEL representative structures whereas CTD was placed at the distance of 20.2Å in average. There is no much deviation in the domain-domain movement of both representative structures. The loop region (206-248) which is believed to be controlling substrate entrance was displaced away from the catalytic cavity in both representative structures. Therefore, based on the free energy minimum cluster population, the cluster1 is selected as a best cluster and representative structure of this cluster used for protein-protein docking studies. Domain-Domain distances and angle of ASM structure were shown in Table 2.  Several studies have revealed that the interaction of PKCδ shows the in uential effect on translocation of ASM from lysosomes to the plasma membrane. In order to understand the molecular level interacting mechanism, the selected representative structures of PKCδ and ASM were used for protein-protein docking. The PKCδ and ASM formed complex with an interface area of 1369.5Å2, constituted by 50 residues of PKCδ and 42 residues of ASM with − 15.4 kcal/mol solvation free energy gained upon interface formation. Both polar and non-polar residues collectively form 16 hydrogen bonds and 5 salt bridges with less hydrophobic contacts. The residues R643, C509, G356, G358, L644, S506, D383, D477, N478, F637 and S645 of PKCδ forms strong hydrogen bonds with the residues S192, Q407, Y500, S501, G502, S503, A176, K186 and R412 of ASM (Fig. 7). The residue S645 of PKCδ which is proven to be an auto phosphorylated residue formed a hydrogen bond with S503 residue of ASM by the contact distance of 3.2Å. This S503 residue is a neighboring residue of S504, which is the site for phosphorylation of ASM and thus these interactions are essential for making phosphorylation process possible.

Structural stability analysis of the bio molecular complex
To check stability of ASM upon complex formation with PKCδ, the complex structure was subjected to molecular dynamics simulation for about 80 ns MD run based on the earlier mentioned molecular dynamics simulation protocol. The backbone RMSD pro le of complex was generated and analyzed. Initially, RMSD uctuation was raised up to 1nm, and was continued up to 35 ns time scale (Fig. 8A). After 35 ns time scale, the RMS uctuation restrained of about 0.2nm uctuation range, which continued till 80 ns with a minimum RMSD uctuation. Additionally, the hydrogen bond pro le shows that 5-10 hydrogen bonds were well maintained in the equilibrated region of 50-80ns MD run describing the stable conformation of the complex (Fig. 8B). Moreover, the FEL analysis showed only one free energy minimum cluster based on the rst two principle components. The representative structure from this free energy minimum cluster is extracted and displayed in Fig. 9, which indirectly claims that the transition of the complex was restricted as the stability is maintained upon complex formation. In addition, the interacting complex was subjected to MMPBSA calculation to nd the binding free energy of protein-protein complex structure in the equilibrated region of MD simulation. Accordingly, best complex shown the binding energy of -585.613 +/-222.974 kJmol − 1 (Table 4).  [55], thus ASM functional activity is very essential. However, its translocation to plasma membrane and subsequent ceramide elevation leads to AD. Hence, we aim to identify inhibitor molecules that could restrict ASM translocation. As an allosteric site, ASM interface that interacts with PKCδ was docked with a library of ligands and the docking results were analyzed. Among the library of compounds, the top four hits were selected based on binding mode, XP score and glide emodel scores. Zinc database code of each ligand and their structures were displayed in Table 5.  Fig. 10B. A recent study has revealed that gentamycin, an antibiotic used for treatment of many anti-microbial infections, targets ASM in cancer cells [57]. This report gives a hope that ZINC71928291 (neomycin), a compound of antibiotic class binding at allosteric site of ASM could be potential lead molecule for the evaluation of its e cacy in AD.  (Fig. 10D). We recognized that ZINC71773625 (rapastinel) is a potent anti-depressant agent, and act as an allosteric modulator of Nmethyl-D-aspartate receptor [58]. Based on our docking results and existing evidence, we believe that this lead molecule could also be a possible allosteric modulator of ASM. For the four complexes, the interacting residues were listed in the Table 6. (rapastinel) showed interactions with residues of ASM, which are essential to make ASM-PKCδ. The existing pro-apoptotic, anti-microbial and anti-depressant activities of ZINC85551993 (3'-Sialyllactose), ZINC71928291 (neomycin) and ZINC71773625 (rapastinel), respectively, and the proven apoptotic role of ASM together supports a possibility that these molecules may act as potential allosteric inhibitors of ASM. Taken together, the four lead molecules may developed as allosteric modulators of ASM for the hope of AD therapeutic implications. and Zinc nger2 domain (cyan). The catalytic domain contains a protein kinase domain (green) and AGC kinase domain (red). ATP (magenta) binds at ATP binding site. The residues that are shown in spherical shape are believed to be involved in auto phosphorylation.   Free energy landscape of PKCδ: FEL map was generated based on rst two principle components whose cosine content value is less than 0.2. It shows three different clusters named as 1 to 3 based on the population size. Three representative structure extracted from each cluster were shown.

Figure 6
Page 25/27 Free Energy Landscape of ASM: The map was derived based on the rst two principle components whose cosine content values are lesser than 0.2. Two representative structures were derived from each cluster to understand domain movements.

Figure 7
The bio-molecular complex of ASM with PKCδ: The best free energy representative structure of ASM (green) and PKCδ (blue) were docked and the complex is displayed here along with the interacting Page 26/27 residues of ASM and PKCδ.  Free energy landscape of complex structure: The most populated cluster was generated and the representative structure conformation was shown.