Targeting Caspase 8: Using Structural and Ligand-Based Approaches to Identify Potential Leads for the Treatment of Multi-Neurodegenerative Diseases

Caspase 8 is a central player in the apoptotic cell death pathway and is also essential for cytokine processing. The critical role of this protease in cell death pathways has generated research interest because its activation has also been linked with neural cell death. Thus, blocking the activity of caspase 8 is considered a potential therapy for neurodegenerative diseases. To extend the repertoire of caspase 8 inhibitors, we employed several computational approaches to identify potential caspase 8 inhibitors. Based on the structural information of reported inhibitors, we designed several individual and consensus pharmacophore models and then screened the ZINC database, which contains 105,480 compounds. Screening generated 5332 candidates, but after applying stringent criteria only two candidate compounds, ZINC19370490 and ZINC04534268, were evaluated by molecular dynamics simulations and subjected to Molecular Mechanics/Poisson Boltzmann Surface Area (MM-PBSA) analysis. These compounds were stable throughout simulations and interacted with targeted protein by forming hydrogen and van der Waal bonds. MM-PBSA analysis showed that these compounds were comparable or better than reported caspase 8 inhibitors. Furthermore, their physical properties were found to be acceptable, and they are non-toxic according to the ADMET online server. We suggest that the inhibitory efficacies of ZINC19370490 and ZINC04534268 be subjected to experimental validation.


Introduction
Neurological disorders (NDs) are usually due to structural and/or functional losses of neurons and eventually neuronal death [1]. Literally, hundreds of neurological/neurodegenerative diseases have been identified, but some like Alzheimer's disease (AD) and Parkinson's disease (PD) are very well known and have devastating impacts in modern society. NDs tend to share some characteristics such

Results and Discussion
Due to the small activity range differences between training set compounds, we used the HipHop protocol for pharmacophore generation. We assumed that the most active ligand(s) in the training set would bind in the same manner with the active site of caspase 8. We assessed the common features essential for binding using the HipHop module in Catalyst software. The six-molecule training set was used to generate pharmacophore models based on common chemical features ( Figure 1).

Figure 1.
Structures of the training set compounds used for pharmacophore generation [26].

Pharmacophore Modeling: 3D Pharmacophore Generation
Ten pharmacophore models (hypotheses) were generated with a rank score ranging from 152.441 to 150.863 (Table 1). The resulting 10 hypotheses contained seven features viz. five hydrogen bond acceptor-lipid (H), one ring aromatic (R), and one hydrophobic (Z) common to all 10 hypotheses. Of these 10 hypotheses, Hypo-1 (Figure 2A) was selected as it mapped all features of the most active molecule in the six-compound training set. Two out of five H functions maps on the five oxygen with different -OH and C=O functions of the most active compound (NP-LHED-AOMK) ( Figure 2B) from the training set, while the R function was mapped to one of the ring aromatic and hydrogen bond donor functions. The NH group from the amino acid linkage supported the remaining single hydrogen donor functionality. One benzene ring from LHED and LTED was mapped for the R function of this pharmacophore model ( Figure 2C).

Pharmacophore Modeling: 3D Pharmacophore Generation
Ten pharmacophore models (hypotheses) were generated with a rank score ranging from 152.441 to 150.863 (Table 1). The resulting 10 hypotheses contained seven features viz. five hydrogen bond acceptor-lipid (H), one ring aromatic (R), and one hydrophobic (Z) common to all 10 hypotheses. Of these 10 hypotheses, Hypo-1 (Figure 2A) was selected as it mapped all features of the most active molecule in the six-compound training set. Two out of five H functions maps on the five oxygen with different -OH and C=O functions of the most active compound (NP-LHED-AOMK) ( Figure 2B) from the training set, while the R function was mapped to one of the ring aromatic and hydrogen bond donor functions. The NH group from the amino acid linkage supported the remaining single hydrogen donor functionality. One benzene ring from LHED and LTED was mapped for the R function of this pharmacophore model ( Figure 2C).  In the direct hit mask, (1) every feature of the training set molecules were mapped; (0) indicates one or more features were not mapped. In the partial hit mask, (0) indicates every feature of the training set molecules were mapped; (1) indicates one or more features were not mapped. Z; Hydrophobic (Z), H; Hydrogen bond acceptor (H), R; Ring aromatic (R).

Pharmacophore-Based Virtual Screening (PBVS)
The validated hypothesis (pharmacophore model) Hypo-1, was used as a query to screen the ZINC database, which contained 105,480 molecules, using the 'Best Flexible Search' option in DS 2.5. A hit list of 5332 compounds matched the query (Hypo1). Further filtration based on fit values of >3.5 and physiochemical screening resulted in 76 hits. These hits were subjected to visual inspection for proper alignments with Hypo-1, and finally, after performing molecular docking, 10 hits were retrieved from the ZINC database ( Table 2). The final selection from this set of 10 hits was made using Moldock and Re-rank scores and resulted in the selection of two compounds (ZINC19370490 and ZINC04534268) for further study.  In the direct hit mask, (1) every feature of the training set molecules were mapped; (0) indicates one or more features were not mapped. In the partial hit mask, (0) indicates every feature of the training set molecules were mapped; (1) indicates one or more features were not mapped. Z; Hydrophobic (Z), H; Hydrogen bond acceptor (H), R; Ring aromatic (R).

Pharmacophore-Based Virtual Screening (PBVS)
The validated hypothesis (pharmacophore model) Hypo-1, was used as a query to screen the ZINC database, which contained 105,480 molecules, using the 'Best Flexible Search' option in DS 2.5. A hit list of 5332 compounds matched the query (Hypo1). Further filtration based on fit values of >3.5 and physiochemical screening resulted in 76 hits. These hits were subjected to visual inspection for proper alignments with Hypo-1, and finally, after performing molecular docking, 10 hits were retrieved from the ZINC database ( Table 2). The final selection from this set of 10 hits was made using Moldock and Re-rank scores and resulted in the selection of two compounds (ZINC19370490 and ZINC04534268) for further study.

Molecular Dynamics Simulation and Molecular Docking
Molecular dynamics simulation of caspase 8 was performed for 20 ns using GROMOS96 53a6 force field. We analyzed the binding efficiencies of ZINC19370490 and ZINC04534268 with caspase 8 for different times (0 ns, 5 ns, 10 ns, 15 ns, and 20 ns). These selected compounds were further tested to confirm their binding affinities. A different scoring function (the gold fitness score) was used to reconfirm their binding affinities with caspase 8. This function scores the binding efficacies of receptor-ligand complexes, and in the present study were used to predict binding affinity [27,28]. The binding efficiencies of both compounds were determined and additional poses were created based on gold fitness scores. The binding efficacies of ZINC19370490 and ZINC04534268 with caspase 8 at different time intervals are shown in Tables 3 and 4.  We found both compounds bound effectively to the active site of caspase 8. The detailed binding mechanisms of the two compounds were further investigated by analyzing binding poses using PyMOL [29]. Figures 3 and 4 depict binding between ZINC19370490 or ZINC04534268 and the active site of caspase 8 at 0 ns, 5 ns, 10 ns, 15 ns, and 20 ns. R179, Y340, R341 and C285 residues in the active site of caspase 8 were found to be most involved in the binding of both compounds. The roles of these key residues have been also described elsewhere [30], and concur with our finding regarding the importance of these residues in the active site. Roles of other important amino acids were also revealed by the present study (Tables 3 and 4). Summarizing, we found caspase 8 tends to form stable complexes with ZINC19370490 and ZINC04534268.
Molecules 2019, 24, x 6 of 11 the importance of these residues in the active site. Roles of other important amino acids were also revealed by the present study (Table 3, Table 4). Summarizing, we found caspase 8 tends to form stable complexes with ZINC19370490 and ZINC04534268.   the importance of these residues in the active site. Roles of other important amino acids were also revealed by the present study (Table 3, Table 4). Summarizing, we found caspase 8 tends to form stable complexes with ZINC19370490 and ZINC04534268.

Thermodynamics Analysis
Thermodynamic analysis produced encouraging results, as binding free energies of ZINC19370490 and ZINC04534268 were either comparable or better than those of previously reported caspase 8 ligands [26] (Table 5). Thermodynamic parameters are important for drug design as they determine the relative binding strengths of protein and non-protein entities. Our MM-PBSA analysis showed the values of different parameters for ZINC19370490 and ZINC04534268 were comparable to those of existing molecules. In addition, we analyzed the solvation free energies of the two proteins to evaluate their water binding affinities, which is a critical consideration when designing molecular therapies. Interestingly, this analysis revealed ligand-binding sites exhibited both positive and negative solvation free energy regions/residues. The number of residues that avoid solvation is in general higher as compared to those that prefer to be soluble. The ligand-binding site of caspase 8 is surrounded by polar residues that pose a barrier to ligand binding.

Physical Properties of Ligands
The physical properties of these ZINC19370490 and ZINC04534268 were found to be compatible with the implementation of pre-clinical studies ( Table 6). Lack of cytochrome P family inhibition is a good indicator of a long circulation half-life, low skin permeability reduces the risk of skin-related injuries, and low gut absorption can be enhanced by techniques such as nanoparticle encapsulation [31]. ZINC19370490 and ZINC04534268 are also non-toxic according to the ADMET online server (https: //preadmet.bmdrc.kr/toxicity/).

Common Feature Pharmacophore Model
Based on our previous studies on approaches to indirect drug design [32,33] and in the field of neurodegenerative disease [33] we decided to focus on caspase 8 inhibitors. Molecules reported in the literature were collated [26], but biological activity data did not exhibit the required 3 log unit variation necessary to develop a quantitative pharmacophore model using the Hypogen module in Catalyst software. Accordingly, we built a common feature pharmacophore 'HipHop' model using six structurally diverse compounds with different inhibitory activities (training set, Figure 1).

Common Feature Pharmacophore Generation
The set of six diverse compounds with maximum activity at the dose of 250µm were used as the training set ( Figure 1). Pharmacophore generation protocol was performed using the HipHop algorithm of Catalyst employed in Discovery Studio 2.0 (DS 2.5) [34]. All training set compounds were drawn/built using ISIS Draw 2.5 and imported into DS 2.5 windows. The CHARMm force field was applied to optimize training set compounds [35]. The conformations of these compounds were generated using the 'diverse conformation generation' protocol of DS 2.5 using default parameters (principal value = 2, maximum omit feature = 0). Accordingly, the most active compound was assigned a score of 1 and moderately and less active compounds were assigned a score of 0. The 'Feature Mapping' protocol was run to detect common features in the training set.

Pharmacophore-Based Virtual Screening (PBVS)
The PBVS approach was used to identify probable inhibitors of caspase 8. The validated model of pharmacophore (Hypo-1) was used as a query to search for compounds in the ZINC database using 'Best flexible search' option in DS 2.5. Resulting hits were screened based on fit values of >3.5 and this was followed by additional screening using physiochemical properties. In addition, these hits were subjected to visual inspection for proper alignments with Hypo-1 and finally subjected to molecular docking, after completing the virtual screening process the 10 most potent hits were retrieved from the ZINC database. Two of these 10 hits were selected based on Moldock score and re-rank scores for further study.

Molecular Docking Studies
Molecular docking was performed using the MolDock module in Molegro Virtual Docker (MVD) software [36]. The scoring function of molecular docking in MolDock is based on piecewise linear potentials (PLPs) [37]. A re-ranking method was applied to the highest ranked poses to increase the accuracy of docking. The search algorithm 'MolDock SE' was applied to this analysis, and population size of 50 and a maximum number of iterations of 1500 were set as parameters. Other parameters were kept as default with the number of runs as 10. Since MVD relies on an evolutionary algorithm, repeated docking runs do not result in precisely the same poses and interactions. To address this intrinsic arbitrariness, three successive runs were performed and the best three poses were used to visualize further interactions.

Molecular Dynamics Simulation
The Gromacs 4.6.7 suite was utilized to perform MD simulations [38] with the 53a6 parameter set of gromos96 force field [39]. A dodecahedron box was selected to immerse the protein and SPC216 water model was used to fill the box [40] that extended to at least 1.2 nm from the edge of the protein. Na + and Cl − counter ions were added to neutralize the system up to a physiological salt concentration of 100mM. After solvation and neutralization steps, system was minimized energetically using the steepest-descent method to remove steric clashes between atoms. In addition, the MD simulations were run under NPT (Constant number of particles, pressure and temperature) conditions with the use of Berendsen's coupling algorithm for maintaining constant temperature and pressure. The LINCS algorithm [41] was used to constrain the lengths of all bonds and water molecules were restrained using the SETTLE algorithm [42]. The particle mesh Ewald (PME) was used to treat long-range electrostatic interaction [43]. A 10 Å cut-off was employed for the van der Waals interactions. The MD simulation time-step was 2 fs and coordinates were saved every 2ps. The PRODRG server was used for the parameterization and generation of topology for each ligand [44].

Thermodynamic Analysis
Protein-ligand binding free energies were calculated to monitor the binding intensities of ligands to caspase 8 using MM-PBSA. For the MM-PBSA calculations, the system was prepared by adding hydrogens, and partial charges on ligands were calculated using the theory utilized in the AM1-BCC model. Calculations were performed using AMBERTOOLS16. Solvation free energies and decomposition into individual amino acids were calculated using ProWave web server (https://www.prowave.org/). The physical properties of ligands were calculated using PreADMET [45] and SwissADME [46].

Conclusions
Common feature-based pharmacophore modelling and structure-based docking study were applied to identify novel caspase 8 inhibitors. A pharmacophore model (Hypo-1) was developed and validated and used as a filtering tool to screen the ZINC database. After a series of screenings followed by visual inspection and molecular docking of selected compounds, the two most potent compounds were subjected to molecular docking with caspase 8 at different time intervals (0 ns, 5 ns, 10 ns, 15 ns, and 20 ns) to obtain deeper insight of their binding efficacies.
The binding free energies of these two compounds were comparable or better than those of previously described caspase 8 inhibitors and satisfied most drug candidate criteria. Using a combination of pharmacophore modeling, virtual screening, molecular docking, and finally molecular dynamics simulations of these two proteins at different time intervals and their binding affinities, we conclude that ZINC19370490 and ZINC04534268 be viewed as potential lead compounds.