Molecular modeling of human neutral sphingomyelinase provides insight into its molecular interactions.

The neutral sphingomyelinase (N-SMase) is considered a major candidate for mediating the stress-induced production of ceramide, and it plays an important role in cell-cycle arrest, apoptosis, inflammation, and eukaryotic stress responses. Recent studies have identified a small region at the very N-terminus of the 55 kDa tumour necrosis factor receptor (TNF-R55), designated the neutral sphingomyelinase activating domain (NSD) that is responsible for the TNF-induced activation of N-SMase. There is no direct association between TNF-R55 NSD and N-SMase; instead, a protein named factor associated with N-SMase activation (FAN) has been reported to couple the TNF-R55 NSD to N-SMase. Since the three-dimensional fold of N-SMase is still unknown, we have modeled the structure using the protein fold recognition and threading method. Moreover, we propose models for the TNF-R55 NSD as well as the FAN protein in order to study the structural basis of N-SMase activation and regulation. Protein-protein interaction studies suggest that FAN is crucially involved in mediating TNF-induced activation of the N-SMase pathway, which in turn regulates mitogenic and proinflammatory responses. Inhibition of N-SMase may lead to reduction of ceramide levels and hence may provide a novel therapeutic strategy for inflammation and autoimmune diseases. Molecular dynamics (MD) simulations were performed to check the stability of the predicted model and protein-protein complex; indeed, stable RMS deviations were obtained throughout the simulation. Furthermore, in silico docking of low molecular mass ligands into the active site of N-SMase suggests that His135, Glu48, Asp177, and Asn179 residues play crucial roles in this interaction. Based on our results, these ligands are proposed to be potent and selective N-SMase inhibitors, which may ultimately prove useful as lead compounds for drug development.

Background

Sphingolipids, mainly sphingomyelin (SM) and glycosphingolipids, are essential constituents of mammalian cells, where they can be found predominantly in the outer leaflet of the plasma membrane [1]. Products of SM metabolism, such as ceramide, sphingosine, sphingosine 1-phosphate, and SM itself, have been recognized as potentially important signaling molecules, which participate in transduction pathways involved in regulation of several processes in mammalian cells. Sphingomyelinases (SMases; EC3.1.4.12) are actually sphingomyelin phosphodiesterases (SMPDs) that catalyze hydrolysis of membrane SM to form ceramide [2]. Ceramide has been suggested to play important roles in cellcycle arrest, apoptosis, inflammation, and eukaryotic stress responses [3]. Production of ceramide by hydrolysis of sphingomyelin activates proline-directed protein kinases, which might be responsible for activation of phospholipase A2 (PLA2) [4, 5]. Production of arachidonic acid by PLA leads to the generation of proinflammatory metabolites [6]. Currently, five types of SMases have been identified and classified according to optimum pH and metal ion dependence ­ the ubiquitous lysosomal acid SMase, the zinc-dependent secreted acid SMase, a neutral Mg2+­dependent SMase, a neutral Mg2+­independent SMase and alkaline SMase [7-9]. Of these, Mg2+­dependent neutral SMase (N-SMase) has emerged as the major candidate for stress-induced ceramide responses. Neutral Mg2+­dependent sphingomyelinases are integral membrane proteins in mammals, and soluble proteins in bacteria. The mammalian sphingomyelinases are thought to play a key role in sphingolipid metabolism and there is increasing evidence implicating SM metabolism in cell signaling, cell proliferation and apoptosis [10-12]. Sequence analysis of these proteins and other eukaryotic N-SMases revealed that they contain either Nor C-terminal extensions encoding predicted membranespanning regions, which will localize the proteins to membranes and promote interfacial catalysis. Subsequent research found that N-SMase activity was induced by a variety of stimuli including cytokines, cellular stresses such as UV light and chemotherapeutic drugs and pathological stimuli like amyloid-β peptides and lipopolysacchride [13-15]. Because of its role in the regulation of important cellular processes, and the continuous availability of substrate SM, the activity of N-SMase must be strictly controlled in cells [16]. Studies addressing the activation of N-SMase by external receptors have strongly focused on the 55 kDa receptor for tumor necrosis factor-alpha (TNF-α) [17, 18]. A region of this TNF-α receptor, the N-SMase activating domain (NSD) adjacent to the death domain, is specifically required for activation of N-SMase through binding of a novel adaptor protein - the factor associated with neutral sphingomyelinase activation (FAN) [19]. The importance of FAN was subsequently confirmed in studies showing that overexpression enhances TNF-stimulated N-SMase activity [20, 21]; therefore, it is likely to play a prominent role in the regulation of major inflammatory cellular responses. In view of the biological importance of N-SMase, we focused our attention to predict the three-dimensional structure of this enzyme and its interacting proteins, specifically TNF-R55 and FAN, using fold recognition and threading methods, with the ultimate aim to gain new insight into their structure and function. This will enable us to understand the catalytic behavior of this enzyme as well as the molecular basis of its activation and regulation.

Methodology

The detail methodologies are provided in supplementary data (available with authors).

Discussion

Model of human N-SMase

In the absence of experimentally determined protein structures, homology-based models may serve as starting points for investigation of sequence-structure-function relationships between diverged enzymes. Although homology-modeled structures may often not be accurate enough to allow characterization of protein-protein or protein-inhibitor interactions at the atomic level, they can suggest which sequence regions or individual amino acids are essential components of the binding surfaces [ http://www.biomedcentral.com/1472-6807/5/2]. However, homology modeling requires a homologous template structure to be identified and the sequence of the protein of interest (the target) to be correctly aligned to the template. The primary protein sequence of human N-SMase was retrieved from the UniProt database (accession number: O60906; Sphingomyelin phosphodiesterase 2, Gene name SMPD2), and the sequence was submitted to NCBI BLAST against the Protein data bank (PDB), in order to identify a template for structure prediction. This analysis revealed the presence of an ENDO/EXO phosphate domain and two transmembrane regions at the Cterminus of N-SMase, which was further confirmed by the SMART server [ http://smart.emblheidelberg.de/help/smart_glossary.shtml]. However, it did not identify any structural templates with sufficient sequence homology and coverage. Thus, in order to identify a template structure for modeling of N-SMase, we used a threading approach, which allows assessing the compatibility of the target sequence with available protein folds based not only on sequence similarity but also on structural considerations (e.g. secondary structure elements). In comparative modeling techniques, the best template is usually chosen based on the percentage of sequence identity between the target and the template proteins. If the overall similarity with the template is low, then the selection can be quite a difficult task as different algorithms may give different hits which may also include false positive results. In this situation, considering results from different servers may provide a more reliable indication. Therefore, the N-SMase sequence was submitted to the GeneSilico metaserver and, additionally, to FUGUE. Out of nine individual servers covered by GeneSilico, eight reported bacterial neutral sphingomyelinase C from Listeria ivanovii (PDBID 1ZWX). Two other templates (2DDS and 2DDR) are also reported by several servers. Despite only 18 % sequence identity with N-SMase, 1ZWX outperformed the rest of the proposed templates; most importantly, the catalytic residues of the two enzymes aligned well with each other. This was not the case for the other templates (2DDS or 2DDR) in which not all of the catalytic residues aligned with catalytic residues of the target sequence. In order to assess the validity of our selected template, we applied the FUGUE program which is able to detect distant homologs by profile matching between the target sequence and a library of known structures. Indeed, FUGUE reported 1ZWX as a template. Therefore, the single consensus template returned by almost all servers (bacterial neutral sphingomyelinase C) was chosen for subsequent calculations. This template gave a reasonably low E-value and, at the same time, provided good coverage of the target sequence (Table 1, see Table 1).

The length of the N-SMase sequence is 423 amino acids, out of which we were able to model the catalytic domain (ENDO/EXO phosphate domain) residues from 1 to 281. Sequence searches and alignments have shown that N-SMases across the phyla share the same basic residues, and domain analysis of the N-SMase family reveals a conserved catalytic domain [10]. The other two transmembrane regions at the Cterminus of N-SMase (residues 282 to 423) were not modeled due to unsatisfactory alignments. Assessment of the reliability of the target-template alignment was achieved by comparing the secondary structure elements of 1ZWX with the consensus predicted secondary structure of the N-SMase catalytic domain. The target-template alignment is shown in Supplementary Figure 1 (available with authors). Fifty models comprising all non-hydrogen atoms were generated using Modeller9v7, structure number 49 was selected as the best model based on the lowest MOF and DOPE score calculated by the program; its loops were optimized using MODLOOP. The alignment shows that the conserved catalytic residues of N-SMase (Glu48, His135, Asp184, and Asp237) aligned well with the catalytic residues of 1ZWX (Glu88, His179, Asp239, Asp282). These include a histidine (His), which is important in acid-base function, a glutamate (Glu) which is essential for Mg2+ binding, and several aspartate (Asp) residues for substrate recognition [10, 11]. It has been reported by Yamada et al. (1988) that the presence of Mg2+ ions interacting with Glu residues increases the activity several fold and also enhances the stability of the protein [22]. Therefore, we docked Mg2+ ion in the modeled structure using Molegro Virtual Docker (MVD) software [ http://www.molegro.com/products.php]. The amino acids involved in metal ion binding of N-SMase are Glu48, His135, His271, and Asn14.

The final N-SMase model, which is illustrated in Supplementary Figure 2 (available with authors), was evaluated using PROCHECK, Verify3D and PROVE. Ramachandran plot generated by PROCHECK showed that 88.3% of the residues were located in the most favored regions, 11.3% in additional allowed regions another 0.4% in generously allowed regions and there were no residues in disallowed regions. In the Verify3D analysis, it was found that 83.3% of the residues scored more than 0.2. For the quality of the model to be considered satisfactory, it is expected to have a Verify3D score greater than 80%. Analysis of the entire structure using the PROVE program gave an average Z-score of 1.783, indicating reasonable packing quality. In PROVE, Z-scores above 4.0 and below ­4.0 suggest the presence of significant errors in the structure in terms of packing of buried atoms. Together, the PROCHECK, Verify3D and PROVE evaluation indicates that the homology model of N-SMase is plausible in terms of prior knowledge, as has to be expected for a valid model. Further, to investigate how well the modeled structure matches the X-ray structure of the template, the N-SMase homology model and the crystal structure of bacterial neutral sphingomyelinase C from Listeria where superimposed with respect to their backbone atoms, yielding and RMS distance of 2.17 Å. We conclude that, despite low overall sequence identity with the template, our N-SMase model meets the quality criteria commonly applied to homology modeled structures, and therefore is considered suitable for further investigation, such as in silico complex formation with physiological binding partners.

Model of TNF-R55 NSD

The domain within the cytoplasmic portion of TNF-R55 responsible for initiating this N-SMase pathway has been functionally mapped to a small region adjacent to the death domain and designated NSD (N-SMase activating domain) [19]. The NSD comprises of 11 amino acids from pos. 338 to 348 (LQKWEDSAHKP) directly preceding the N terminus of the death domain. Considering the importance of this interaction, we predicted the structure of FAN binding NSD motif using the I-TASSER method. The resulting model is in good agreement with its predicted secondary structure, and it has about 90.6% of amino acid residues in the favored region and about 9.4% of residues in the additional allowed region of the Ramachandran plot. The structure was used for identifying the binding sites for proposed protein-protein interactions. The modeled structure is shown in Supplementary Figure 3 (available with authors).

Model of FAN:

The cytoplasmic domain of TNF-R55 is responsible for TNFinduced activation of N-SMase but there is no direct association between TNF-R55 and N-SMase, which could provide hints on the activation mechanism. Indeed, an adaptor protein termed FAN (factor associated with N-SMase activation) has been reported to couple the TNF-R55 to N-SMase [19]. FAN represents a member of the WD40 (or WD-repeat) family that includes mainly regulatory proteins. The WD-repeat motif consists of a conserved core of 23-41 amino acids, usually bounded by Gly-His (GH) and Trp-Asp (WD); several of these motifs arranged in tandem typically give rise to so-called β- propeller folds which are mainly implicated in protein-protein interactions. Since FAN appears to be crucially involved in TNF-induced activation of the N-SMase pathway, we predicted its molecular structure by homology modeling, with the aim to investigate complex formation with relevant partners.

The amino acid sequence of human FAN was retrieved from the UniProt database (Primary accession number: Q92636, Factor associated with neutral sphingomyelinase activation, Gene name FAN). Sequence analysis indicates that the protein contains three domains, namely a BEACH domain, a GRAM domain and the C-terminal WD repeat domain; the latter has been shown to be responsible for N-SMase activation [19]. Therefore, we decided to focus on WD repeat domain residues 703-916. The NCBI BLAST server was used to find homologous proteins for comparative modeling. Based on these results, WD repeat-containing protein 5 from Escherichia coli (PDB code: 3EMH) was selected as a template, with 27.8% identity and 48% similarity between sequences; for the alignment please refer to Supplementary Figure 4 (available with author). The final model (Supplementary Figure 5 available with authors) is in good agreement with its secondary structure prediction, and it has about 85.9% of amino acid residues in the favored region and about 13% of residues in the additional allowed region of the Ramachandran plot. The model was further evaluated with Verify3D and it was found that 85.7% of the residues scored more than 0.2; analysis of the entire structure using the PROVE program gave a Z-score of 1.584.

Molecular dynamics simulation

To gain insight into stability and dynamics of the predicted structures, explicit solvent MD simulations were performed for all three proteins (N-SMase, FAN, and NSD) for 3 ns. The total energy of the system and the RMSD from the starting structure are essential to evaluate the sustainability and convergence of any MD simulation. For the three proteins considered here, Supplementary Figure 6 (available with authors) presents the RMSD of Cα backbone atoms with respect to the initial structure. The graphs clearly show that the RMSD mostly did not exceed 3.5 Å for N-SMase and FAN (panels 1A and 2A, respectively), and about 3.0 Å in case of NSD (panel 3A). These results suggest that a relatively stable conformation of all three proteins is achieved through MD simulation. Furthermore, from the analysis of the MD trajectories we note a gradual decrease in the potential energy of the models (panels 1B, 2B, 3B), which indicates that the structures are energetically stable during MD.

In order to identify major motion of the protein Principle Component Analysis (PCA) tool is used. The method is based on the calculation and diagonalization of the covariance matrix of atomic fluctuation, which yields a set of eigenvalues and eigenvectors. The covariance matrix of N-SMase, FAN, and NSD were constructed based on a trajectory of 1501 frames with dimension 2529 for N-SMase, 2682 for FAN and 99 for NSD. These dimensions were given by 843 backbone elements for NSMase, 894 for FAN, and 33 for NSD, with eigenvalues sums of 10.7083, 11.1347 and 1.70545 nm2. These eigenvalues were plotted against the eigenvectors of each protein yielded steep eigenvalue curves (Supplementary Figure 7 available with authors). Thus the largest eigenvalues describe the dominant motion within the simulation and it reveals that 90% of the backbone motion is covered by the first 20 eigenvectors. These result shows that most of the internal motion of the protein is restricted to a subspace with very small dimension. Furthermore, the motion along the principle modes was analyzed by projecting the atomic trajectories onto the corresponding eigenvectors. Supplementary Figure 8 (available with authors) panels 1A, 2A, 3A) shows the first three, the tenth and the twentieth projection of the backbone trajectory on the eigenvectors obtained from backbone covariance matrix are plotted against the time for each of N-SMase, FAN, and NSD. It is clear from the figure that the motion of the backbone reached their equilibrium fluctuations in the first ten eigenvectors. Supplementary Figure 8 (available with authors) (panels 1B, 2B and 3B) shows the trajectory projected on the planes defined by two eigenvectors (the tenth and twentieth) from the backbone coordinate covariance matrix of N-SMase, FAN, and NSD are strongly correlated and belong to same set. This indicates that the trajectories fill the expected range almost completely and there is no high projection observed far from the diagonal. The result obtained shows that modeled structures show prolonged stability in the time scale of 3 ns MD simulation, which in turn favors the selection of these structures for further analysis.

Protein-protein docking

Protein-protein interaction is critical to understand many biological processes such as the regulation of protein activities. It has been reported that human N-SMase is activated by TNF through FAN [19]. To gain structural insight into the role of different residues in this interaction, we performed proteinprotein docking for TNF-R55 NSD, FAN and N-SMase; in doing so, we made use of available experimental evidence for a rough pre-orientation of the proposed binding surfaces. We used different docking servers (i.e. PATCHDOCK, FIREDOCK, GRAMMX, CLUSPRO) to obtain possible complexes. It has been reported that TNF-R55 NSD binding to FAN leads to a conformational change in the FAN protein, enabling it to mediate activation of N-SMase. Therefore, we first docked TNFR55 NSD with FAN and, using this assembly, we investigated association with N-SMase to obtain the ternary complex. The best cluster returned by the PATCHDOCK server was chosen based on its calculated binding energy (­62.72 kcal/mol), and was further refined by FIREDOCK.

For this protein complex (Figure 1), we evaluated the interacting residues, the attractive and repulsive van der Waals (VdW) energies, and atomic contact energies (ACE, Table 2, see Table 2). These predicted functional sites can be considered one of the most essential structural regions to increase the enzyme activity by conformational changes without compromising the integrity of the TNF-R55 NSD-FANN- SMase complex. Analysis of the protein complex using Dimplot (not shown) confirms that FAN couples TNF-R55 and N-SMase; this is in agreement with experimental data, suggesting a crucial role in mediating TNF-induced activation of the N-SMase pathway, which in turn regulates important mitogenic and proinflammatory responses.

Figure 1

Model of the TNF-R55 NSD-FAN-N-SMase ternary complex. The molecules are displayed in cartoon sketch and the interaction surface is in mesh representation. The figure was prepared using Discovery Studio 2.0.

Molecular dynamics of the TNF-R55-FAN-N-SMase complex

The TNF-R55-FAN-N-SMase complex has been subjected to molecular dynamics using the same protocols as described in the Methods section, but for a total duration of 1000 ps. Trajectory analysis confirms the stability of the complex. Specifically, protein conformations in the predicted interaction site were found to be stable during the simulation with respect to starting structures. Conformational changes during the simulation were monitored by Cα RMSD computed for the complex (Supplementary Figure 9 available with authors) Furthermore, the gradual decrease in potential energy of the model from ­165,000 kJ/mol to ­166,000 kJ/mol indicates that the complex is energetically stable during MD.

Protein-ligand docking studies

To study the molecular basis of interaction and affinity of glutathione (GSH), which inhibits the N-SMase in vitro [54], a library of selected ligands were docked into the active site of the enzyme. The docking results for these ligands are given in Supplementary Table 3. Out of 15 candidate complexes generated by XP docking, potential hits were selected by visual inspection and studied in detail. The two-dimensional molecular structures of the ligands are displayed in Figure 2. Docking of Ligand 900 into N-SMase suggested hydrogen-bond interactions with four residues, Glu48, His135, Asp177 and Asp179 (Supplementary Figure 10a available with authors). This implies a competitive mechanism of inhibition as these residues are part of the conserved catalytic site. The various residues involved in van der Waals interaction are His181, Glu139, Ala138, His 135, Trp16, which provide stability to the complex. The observed Glide docking and emodel scores for ligand 900 were ­9.90 and ­53.30, respectively. Six residues of the N-SMase binding site, His271, Ser209, Gln224, Asn179, Asp177, and Glu39, establish hydrogen bonding interaction with Ligand 901, in addition to the participation of Val234, Pro231, Lys228, Leu227, and Ser209 in van der Waals contacts (Supplementary Figure 10b). This compound mainly binds to the substrate recognition site and displays a different binding mode with Glide docking and emodel scores ­9.26 and ­68.26. The Glide emodel value of Ligand 902 is very low at ­71.73 (Table 3, see Table 3), which suggests an energetically favorable interaction with the protein. Ligand 902 suggests its binding to five residues of N-SMase, forming five hydrogen bonds (Supplementary Figure 10c available with authors). The lengths of these bonds are less when compared to other ligands, which supports the concept of a stronger interaction. Among the residues involved, Glu48, Asn179 and Asp177 play important roles in catalysis and substrate recognition. The docking of Ligands 903 and 904 suggests interaction with five and three residues of N-SMase, respectively (Supplementary Figure 10d and 10e available with authors). Ligand 903 binds to His271, Tyr140, Ser113, Asp110 and Glu48, forming a total of five hydrogen bonds; in case of Ligand 904, three bonds are observed. Nevertheless, the emodel score of Ligand 904 is lower than for Ligand 903 (­68.15 vs.­ 57.68). The well-known inhibitor GSH, i.e. Ligand 905, was docked in the active site of N-SMase. The docking results show that it forms four hydrogen bonds with the N-SMase binding site (Supplementary Figure 9f available with authors). The residues involved are Arg151, Tyr140, Glu139 and Hie271; Ile236, Ser209, His181, Glu139, Ala138, His135, Lys115 and Trp50, which are engaged in van der Waals contacts, thus enhancing the stability of the complex. The observed Glide docking score is ­5.72 and emodel score is ­47.92. As a control we re-docked co-crystal structure of bacterial neutral Sphingomyelinase (PDBID 2DDT) using Glide docking tool. The co-crystal ligand 2-(N-Morpholino)-Ethanesulfonic Acid is docked into the 2DDT structure. The supplementary Figure 11 (available with authors) shows that Glide docking tool predicted the correct binding pose of the ligand where it has three hydrogen bond interactions exactly as in the experimental results. The Glide docking score and emodel score are ­6.416 and ­43.620 respectively. This docking result shows evidence that Glide can predict ligand conformations quite well. The ligands predicted in this study can be suitable inhibitors for NSMase.

Figure 2

Two dimensional molecular structures of the ligands identified in our docking studies with N-SMase, along with their IUPAC names. The figure was prepared using Marvin View 5.1.

Conclusion

The main objective of our work was to construct a model of NSMase, as knowledge of the three-dimensional structure is essential for a better understanding of the functional mechanism. After comparative modeling, structural refinement by 3 ns MD simulation improved the general quality of the structure as assessed by various model evaluation methods. Despite the low sequence identity between human N-SMase and its homolog from L. ivanovii, which was used as the template, the final model therefore does have acceptable profiles in our structural analysis. Furthermore, we have predicted the structure of a ternary complex which provides a better understanding of the molecular interactions with FAN and TNF-R55 NSD. We hope that our models will inspire new experimental efforts in this area; specifically, the critical residues identified by our in silico approach can be readily tested for their biological relevance. From the analysis of ligand docking, we conclude that the N-SMase residues most crucial for interaction with these compounds are His135, Glu48, Asp177, and Asn179. We believe that the ligands reported here may emerge as potent inhibitors of N-SMase, possibly opening up a new avenue towards management of excessive inflammatory responses in humans.