Screening of FDA approved drugs for finding potential inhibitors against Granzyme B as a potent drug-repurposing target

Journal Pre-proof Screening of FDA approved drugs for finding potential inhibitor against Granzyme B as a potent drug-repurposing target Saima Ikram, Jamshaid Ahmad, Serdar Durdagi PII:

S1093-3263(19)30493-0

DOI:

https://doi.org/10.1016/j.jmgm.2019.107462

Reference:

JMG 107462

To appear in:

Journal of Molecular Graphics and Modelling

Received Date: 29 June 2019 Revised Date:

12 September 2019

Accepted Date: 29 September 2019

Please cite this article as: S. Ikram, J. Ahmad, S. Durdagi, Screening of FDA approved drugs for finding potential inhibitor against Granzyme B as a potent drug-repurposing target, Journal of Molecular Graphics and Modelling (2019), doi: https://doi.org/10.1016/j.jmgm.2019.107462. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Screening of FDA approved drugs for finding potential inhibitor against Granzyme B as a potent drug-repurposing target Saima Ikrama,b, Jamshaid Ahmada,*, Serdar Durdagib,*

Screening of FDA approved drugs for finding potential inhibitor against Granzyme B as a potent drug-repurposing target Saima Ikrama,b, Jamshaid Ahmada,*, Serdar Durdagib,* a b

Center of Biotechonology & Microbiology, University of Peshawar, Pakistan

Computational Biology and Molecular Simulations Laboratory, Department of Biophysics, School of Medicine, Bahcesehir University, Istanbul, Turkey

Abstract Granzyme B is one of the best-characterized and extensively studied members of cytotoxic lymphocytes (CL) proteases. Initially is thought to be involved in eliminating virally infected or cancerous cells by using a specialized mechanism through which they are internalized into target cells. In the last decade, however this dimension has changed as there are several reports show that not only CL but also other immune cells can also synthesize Granzyme B. This leads to the possibility of the presence of these proteases in extracellular environment. Being active protease, it then raises the possibility of damaging host tissues as evident from the available reported literature. In many instances, Granzyme B is directly involved in pathogenicity, however in others, it contributes to the disease severity as their over expression makes the clinical situation quite worse which ultimately leads to the chronic state of the disease. Serine protease inhibitor-9 is a natural known intracellular inhibitor of Granzyme B, however there is less data available about the potential inhibitors that can regulate their activity in an extracellular environment. Current study is an effort to identify potential novel inhibitors of Granzyme B. For this aim, drug repurposing study was performed. Around 7900 FDA approved drugs were screened using both ligand- and target-driven approaches. Initially all molecules were docked using induced fit docking (IFD) approach and selected 318 high-docking scored molecules were used in initially short (1-ns) molecular dynamics (MD) simulations. Based on MM/GBSA binding free energy calculations, 6 compounds were selected and used in long (100-ns) MD simulations. These compounds were then used in binary QSAR analysis. Therapeutic activity potentials of studied compounds were investigated by Clarivate Analytics’s MetaCore/MetaDrug platform which uses binary QSAR models. It

1

is developed based on manually curated database of molecular interactions, molecular pathways, gene-disease associations, chemical metabolism and toxicity information. Results of selected compounds were compared with a positive control molecule. Current drug repurposing study is a step ahead in finding a potential lead compound by screening database of FDA approved molecules. We have identified inhibitors (Tannic acid, Mupirocin, Phytonadiol sodium diphosphate, Cefpiramide, Xenazoic acid) with a potential of obstructing the active site of Granzyme B.

Key words: Cytotoxic lymphocytes, Granzyme B, serine proteases, GASPIDs, molecular docking, MD simulations, binary QSAR, drug repurposing

*E-mails: [email protected] (JA); [email protected] (SD)

2

Introduction Immune system provides protection against diseases by identifying and killing pathogens, infected and tumorous cells. Variety of immune cells like; cytotoxic lymphocytes (CL) including cytotoxic T lymphocytes (CTL) and natural killer (NK) cells, mast cells and neutrophils, mediate this defense system. These immune cells express a variety of serine proteases, known as granule associated serine peptidases of immune defense (GASPIDs) [1, 2]. Best-characterized GASPIDs include; neutrophil elastase, cathepsin G, proteinase3, granzyme A, B & H, mast cell chymase and tryptase. However, little is known about other members. Although GASPIDs are highly related, they may show different substrate specificities implying different biological functions. GASPIDs which are synthesized by CL are known as granzymes. In human, there are five granzyme genes (Gzm-A, B, H, K and M). Based on their specificities these five human granzymes can be grouped into: tryptases (cleavage after basic residues i.e., GzmA and K [3, 4]: asp-ases (cleavage after acidic residues i.e., GzmB) [5]; and chymases (cleavage after hydrophobic residues i.e., GzmM and H) [6-8]. They are located on three chromosomal loci in human genome [1, 9]. Granzymes are synthesized as single domain serine proteases of approximately 240 to 250 amino acids length, encompassing a leader sequence used for trafficking to ER, an N-terminal pro-dipeptide, followed by a mature protein start sequence i.e. I-I-G-G and catalytic triad that are involved in catalytic activity of the granzymes i.e., His57, Asp102, Ser195 (chymotrypsin numbering). At the sequence level, granzymes contain a family specific conserved sequence P-H-S-R-P-Y-M-A, starting at 9th to 16th amino acids after the mature N-terminus [1]. These granzymes perform different intra– and extracellular functions however some of the key roles are described. Three core intracellular roles, dependent on perforinmediated delivery into target cells, have been described: Cytotoxicity, accessory functions and cytokine signaling (Figure 1). Cytotoxicity is the direct killing of cells by activating intrinsic apoptotic pathway and it is mostly induced potently induced by human GzmB, which cleaves the pro-apoptotic molecule Bid (to activate the mitochondrial pathway of apoptosis) as well as directly activating caspases [10-12]. Accessory functions are best described by GzmH, which degrades a virally encoded

3

GzmB inhibitor, thus enabling GzmB to function [13]. Cytokine signaling is carried out by GzmA, which is involved in a poorly characterized pathway that results in the maturation and release of immuno-modulating cytokines [14-16]. Extracellular functions, independent of perforin, have been reported for GzmA and B and involve cleavage of extracellular matrix and cell surface receptors [17], possibly to facilitate CL migration (Figure 1). The roles of GzmK and M are currently unclear, however it has been shown to play a role in generating inflammatory response [18]. In addition, these extracellular proteases also activate cytokines.

Figure 1: Different roles performed by granzymes. Granzymes are internalized into a target cell thorough granule exocytosis method. Cytotoxic role is shown when Granzyme B is entered into a target cell through immunologic synapse by the pore forming protein i.e. perforin. The pathogens are also phagocytosed by immune cells and when fused with lysosomes containing granzymes, ultimately resulting in phagolysosome, thus killing the pathogens (2). Granzymes are also involved in the activation of pro-cytokines to generate

4

inflammatory response (3) and cleavage of extracellular matrix to facilitate the moment of immune cells (4). Regulation of Granzymes As granzymes are the key mediators of various functionalities associated with CL, it is important to regulate their activities. Otherwise unregulated activities of granzymes may result in various complications. Different mechanisms are implied for this regulation however one of the mechanisms is to regulate their activity at post-translation level by a family of serine protease inhibitors known as serpins. The interaction results in the formation of reversible protease-inhibitor complexes. This facilitates the docking of reactive center loop (RCL) that acts as a pseudo substrate, into the active site of the protease. This either results in the inhibitory complex thereby inactivating the protease or it might follow the non-inhibitory pathway, in which serpins is cleaved into in its inactive form with the release of a protease [19]. Among many serpins, only three have been reported to function as granzymes inhibitors [19]. SERPINB9 (PI-9), a member of clade B serpin, is one the best-characterized serpin. It is a potent inhibitor of human GzmB [19, 20]. The search for finding inhibitors is underway however it is in preliminary stage and further investigation needs to be performed for the determination of novel inhibitors, which can then be possibly, used for the regulation of GzmB. One of its inhibitor which is tetrapeptide aldehyde was docked to the binding pocket of the target and the crystal structure of the complex was also resolved. It confirmed the preference for acidic residues as the preferred cleavage at P1 site [21]. Extracellular Granzyme B GzmB is one of the best-characterized members of granzymes family and it is a general consensus that this protease is involved in cytotoxicity besides performing many other roles in immune system. Although the entry of GzmB into a target cell is precisely regulated, in certain instances, GzmB escapes the immunological synapse which then become a part of extracellular milieu. The granule exocytosis pathway is a refined mechanism through which this protease internalizes into a target cell that involves the use of perforin which is a pore forming protein, co-synthesized and packed inside azurophilic granule by CL (Figure 1). After entry into the target cell, it then performs various intracellular functions. During this internalization, it is further reported that sometime

5

about one third of the GzmB escapes into the outside environment from immunological synapse [22] (Figure 1). It is further observed that number of other immune cells like includes macrophages, CD8+ & CD4+ cells, neutrophils, basophils; dendritic cells also synthesize GzmB [23-30]. It is further noted that no extracellular inhibitor is reported for GzmB, although its shows activity, making proteolytic cleavage of various tissues. Implications of Extracellular Granzymes B in disease conditions Initially it is thought as immune facilitator that these proteases act in a substrate specific manner. However, in recent years the paradigm has now been changed and various studies have revealed over expression of these proteases in autoimmune, chronic inflammation, cardiovascular, graft versus host reactions and other metabolic disorders [31-34]. The over expression of these proteases worsens disease conditions by creating highly proteolytic environment [34]. Myasthenia Gravis (MG), which is one the beststudied neuromuscular autoimmune disease. The antibodies are generated against epitops on the acetyl choline receptors (AChRs) which triggers antibody-mediated autoimmune attack. It is generally observed that abnormalities in thymus is involved in this disease. GzmB is reported to be present in elevated amount in the thymus of MG patients while it is absent in normal individuals. The acetyl choline receptors are reported to be cleaved by GzmB, thus generating neu-antigens against which antibodies generated immune attach is initiated, thereby leading to autoimmune conditions [35].

Similarly, GzmB is also

reported to up-regulate inflammatory response by cleaving cytokines like IL-1α, IL-18. IL-18 is present in normal amount in healthy tissues, however in various skin disease, its elevated level is attributed to the cleavage by excess amount of GzmB in the skin environment, therefore it contributes to the chronic inflammation in various skin disorders like cutaneous lupus erythematosus, psoriasis and dermatitis [36, 37]. GzmB is also involved in the activation of protease-activated receptors (PAR)-1 on the surface of neural cells. This activation then initiates voltage gated channel which in turn results in activation of Notch-1 which results in neuro-degeneration or neurotoxicity, ultimately leading to neuro-inflammatory brain disorder [38]. The increased expression of GzmB is found in heart diseases. The over expression of GzmB is observed in angina pectoris, ruptured coronary plaques. The excess of this enzyme is attributed to thrombolysis in myocardial infarction [39]. As granzyme B is also reported to have a role in cleaving

6

extracellular matrix proteins like laminin, vitronectin, and fibronectin [40], it therefore degrades extracellular matrix proteins associated with heart tissues resulting in atherosclerotic plaque vulnerability [41]. Similarly it is reported to be a major player in various renal and other coronary artery diseases [39, 42]. Based on the reported data about the overexpression of granzymes in various diseased conditions, it then becomes imperative to investigate novel inhibitors, which can possibly be used to inhibit the activity of GzmB in various clinical situations, thereby leading into a potential drugable target. The present study will investigate updated FDA approved drugs database for the identification of novel inhibitors. Such understanding can be of paramount importance in designing and developing novel therapeutic inhibitors with low risk for immunogenicity in humans. Methods The overall methodology used in this study was represented and summarized at the Figure 2.

Figure 2: The workflow diagram shows the steps performed in the current study. Protein structure was downloaded from RCSB PDB database and ligands are downloaded from NIH Chemical Genomics Center (NCGC) Pharmaceutical Collection (NPC) database. Proteins and ligand are then prepared for the docking studies. After docking studies, molecular dynamics (MD) simulations and post-processing MD analyses were then

7

performed. Selected molecules were then passed through MetaCore / MetaDrug software to check the therapeutic activity profiles.

Protein Preparation

Crystal structure of human GzmB (PDB ID, 1FQ3) was taken from protein data bank (PDB) and prepared using the protein preparation module of Maestro molecular modeling package. Hydrogen atoms were added, bond orders were assigned and protonation states of residues was determined at the physiological pH 7.4 [43, 44]. The prepared protein structure was then used in geometry optimization using OPLS2005 force field [45]. Ligand Preparation

Around 7900 FDA approved drugs database was downloaded from the NIH Chemical Genomics Center (NCGC) Pharmaceutical Collection (NPC) database. The downloaded database was prepared using the LigPrep module of the Maestro molecular modeling package with the OPLS2005 force field [45]. After the ligand preparation, total number of compounds were reached to ~40000 due to different combinations of enantiomers and tautomers of prepared molecules. Epik was used to apply the protonation states of the ligands at pH 7.4 [46]. Selection of Binding site

GzmB is a serine protease and the active site comprises of catalytic triad that constitutes three key amino acid residues (His57, Asp102 and Ser195, Chymotrypsin numbering) [47]. In addition to catalytic triad there are other specificity conferring residues, located at positions 189, 216, and 226 [47]. GzmB being an asp-ase has a preference for cleavage after acidic residues in a polypeptide chain. Molecular Docking

In the docking studies, various docking algorithms including Glide/SP, Glide/XP and Glide/IFD approaches were used in Maestro with flexible and rigid ligand sampling [48, 49]. A receptor grid box was generated and side chains of certain amino acid residues (i.e., Thr, Tyr) were allowed to rotate to add flexibility to the used target structure. Both flexible and rigid docking algorithms were used to generate the docking poses with realistic docking scores and interaction with catalytic and other active site residues. Binding pocket was created around crucial residues of GzmB protein (His57, Asp102,

8

Ser195). For the generation of different binding poses of ligand, OPLS-2005 force field was used [45]. IFD method provides a better flexibility of both ligand and binding pocket residues throughout the docking. IFD procedure was established in three consequence stages which involves: (i) docking of all studied compounds at the binding pocket of the target structure; (ii) structural refinement of amino acid residues within 5 Å of docked poses (i.e., energy minimization of binding pocket which includes ligand and residues within 5 Å from the docked pose were used in geometry optimization), (iii) refined geometry of the target protein in the previous step were used in re-docking of ligands. While in initial docking, Glide/SP approach was used, in re-docking Glide/XP protocol was performed. Molecular Dynamics (MD) simulations

To analyze the structural and dynamical profiles of ligands at the binding pocket of targeted protein, MD simulations were performed using the Desmond program with the OPLS2005 force field and RESPA integrator. [50] The complex generated from topdocking pose of IFD, was immersed in orthorhombic water box (minimum distance between the box of edges and complex set to 10 Å). Explicit water molecules (SPC solvent model) were used in the preparation of the system and 0.15 M NaCl salt were added to the box in order to neutralize the used systems. In all MD simulations, NPT ensemble at 310 K with Nose-Hoover temperature coupling [51] and at constant pressure of 1.01 bar via Martyna-Tobias-Klein pressure coupling [52] were used. There were no constraints on the generated systems and the initial velocity values are used as default. 2.0 fs time step was used in simulations. In MD simulations, two different production run times were considered: (i) 1 ns (short MD simulations); (ii) 100 ns (long MD simulations). A threshold docking score (-10.00 kcal/mol) was considered for short MD simulations and 318 top-docking scored complexes which have < -10.00 kcal/mol topdocking scores were selected. Since the positive control molecule has -7.35 kcal/mol docking score value, -10.00 kcal/mol as cut-off value may provide safer area for the selection of potent therapeutic compounds. The simulation boxes were prepared individually for these 318 protein-ligand complexes with an in-house script. Molecular Mechanics / Generalized Born Surface Area (MM/GBSA) calculations were performed for the derived trajectory frames throughout the 1-ns MD simulations (100 trajectory

9

frames were used). Based on average MM/GBSA scores 6 compounds were selected for the longer (100-ns) MD simulations. For the MD simulations of selected 6 compounds as well as for the reference compound, simulations are replicated. Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) Calculations.

Schrodinger's Prime module was used to execute MM/GBSA binding free energies calculations of selected ligand-protein complexes throughout the MD simulations. 100 trajectory frames for short MD simulations and 200 trajectory frames for long MD simulations were chosen, respectively. For the prediction of free binding energies of complexes, OPLS2005 forcefield and VSGB 2.0 solvation model were applied [53]. MetaCore / MetaDrug Analysis

Top-docking scored FDA approved drugs which show chemical interactions with catalytic triad residues were selected for the therapeutic activity prediction protocols. These compounds were used in MetaCore/MetaDrug platform. [54] MetaCore/MetaDrug is an inclusive tool for the analysis of compounds and their biochemical and pharmacological behaviors in the field of drug design and development. MetaDrug uses the property of Tanimoto Prioritization (TP) to find the similarity between analyzed compounds and compound sets in the Quantitative Structure Activity Relationships (QSAR) models based on elements found in the structure. These models were prepared with a diverse set of compounds based on experimental evidence of their activity/function on a certain protein of interest, and then tested with validation sets. The accuracy of each model depends on the number of compounds used to create it and can be estimated by the correlation coefficient (R2) and root mean squared error (RMSE), where a higher R2 and low RMSE indicate higher model accuracy. The QSAR model with the highest specificity, sensitivity, accuracy and the Matthews Correlation Coefficient (MCC) was selected in MetaCore/MetaDrug for each tested particular therapeutic activity. The prediction of a therapeutic activity is calculated based on the ChemTree ability to correlate structural descriptors to that property using the recursive-partitioning algorithm. The ChemTree parameters that gave the best results were as follows: path length 5, max segments 3, p-value threshold Bonnferoni 0.99, p-value multiway split 0.99 and number of random trees 50. The training set used in MetaCore/MetaDrug includes molecules that possess the property (positives) and chemicals that do not have such property (negatives)

10

in approximately equal numbers. The marketed drugs were used if their number was greater than 100 in the disease QSAR models; otherwise drug candidates in clinical trials and preclinical compounds with in vivo activity have been added to the training set. Results and Discussion Crystal structure of the human GzmB was used from protein data bank (PDB, 1FQ3) [47]. This protein was expressed in a baculovirus system, crystalized without any inhibitor with 3.1 Å resolution. Binding pocket residues were denoted on the basis of chymotrypsin numbering (catalytic triad residues were His57, Asp102 and Ser195). Other catalytic elements that constitute the oxyanion hole are arranged at the center of the cleft as in trypsin and chymotrypsin. Besides catalytic triad, the other specificity conferring residues are located at positions Thr189, Gly216 and Arg226 that are involved in making specificity pocket [47] (Figure 3).

Figure 3: The crystal structure of human Granzyme B. While red colored region shows the catalytic triad residues (His57, Asp102 and Ser195), blue colored region shows the pocket specificity residues (Thr189, Gly216 and Arg226).

Structural and experimental studies have shown that there are at least six regions of GzmB that are involved in the inhibitor-proteases interactions. The S1, S2 and S3

11

constitute the residues that are involved in catalytic activity, however residues at positions 40, 41, 171, 189, 191, 192, 216 and 226 are also crucial for interactions [47] (Figure 4).

Figure 4: Substrate-protease complex showing specificity pockets in active site of the protease. The specificity pocket is represented by S1 and S1’, corresponding substrate residues that fit into these pockets are represented by P1 and P2’. Before the docking, used human GzmB was prepared using protein preparation module of the Maestro molecular modeling package. Missing side chains were fixed and the protonation states of amino acids at pH 7.4 were generated using PROPKA to mimic physiological conditions [43]. The OPLS2005 force field was finally used for the structural optimization of the protein. Glide/IFD method was used for the estimation of binding free energies as well as binding poses of used FDA approved database which includes around 8000 small molecules. Docking poses of all compounds were analyzed and high docking-scored hits that show chemical interactions with the triad catalytic residues were considered for further analysis. 318 top-docking scored (i.e., docking scores <-10.0 kcal/mole) compounds were selected and used in short MD simulations. Simulation boxes of top-docking poses of 318 complexes are prepared and 1-ns MD simulations were performed using Desmond. Based on MM/GBSA scores (Figure 5), 6 top-MM/GBSA scored molecules were selected for long MD simulations.

12

Figure 5. The graphs show the normalized distribution curve for the measured MM/GBSA scores (kcal/mol) of studied compounds and Z-Scores. Glide/IFD docking scores of these 6 compounds are shown in Table 1. Table 1 also includes docking score of a well-known inhibitor of human GzmB. As it is shown, all selected compounds have better interaction energies compared to positive control. Since reference compound and selected molecules may not have similar size (molecular weight), we also investigated ligand efficiency scores. Ligand efficiency scores are the docking scores per non-hydrogen atoms. Ligand efficiency scores also show that these 6 selected compounds have better interaction energies compared to a known GzmB inhibitor (i.e., reference molecule [55]). (Table 1)

13

Table 1: Docking score of Glide/IFD and ligand efficiency is shown for the 6 compounds which showed high MM/GBSA scores. Table also includes corresponding result for a known GzmB inhibitor (reference molecule [55]). Molecules

2D Structure of molecules

Glide/SP

Glide/IFD

IFD score

(kcal/mol)

(kcal/mol)

per heavy atoms

(kcal/mol) Tannic acid

No pose

-17.54

-0.38

Mupirocin

-4.60

-13.73

-0.39

Phytonadiol

-7.04

-12.45

-0.30

-5.27

-11.14

-0.26

sodium diphosphate

Cefpiramide

14

Xenazoic acid

-5.49

-10.75

-0.38

Vidarabine

-5.98

-10.25

-0.44

Reference

-4.97

-7.35

-0.16

molecule [55]

Long (100-ns) MD simulations were performed for the selected 6 compounds and in order to compare the results, same protocol was also used for the reference molecule [55]. The Root mean square deviation (RMSD) is used to measure the average change in displacement of a selected atoms for a particular frame with respect to a reference frame. For all RMSD and root mean square fluctuations (RMSF) trajectory analysis, 2000 trajectory frames were considered throughout 100-ns MD simulations. RMSD is measured using Cα atoms of the complex structures throughout the MD simulations (Figure S1). As it can be seen from RMSD-time graphs, maximum RMSD value was less than 3.0 Å and all structures were in structural equilibrium after 10-ns. While reference molecule showed fluctuations between 1.5 to 2.2 Å of Cα RMSD values at the first half of the simulations, it is decreased to around 1.5 Å throughout the second half of the simulations. There were slight fluctuations evident from an increase of RMSDs from 2.1 to 2.5 Å over the 100-ns run time in tannic acid, mupirocin and cefpiramide. While average Cα RMSD of protein is found as 1.92 Å when reference molecule was bound, corresponding values were 1.75, 1.66, 1.49, 2.03, 1.91 and 2.16 Å for xenazoic acid, phytonodiol sodium diphosphate, vidarabine, tannic acid, cefpiramide and mupirocin,

15

respectively. Thus, it is found that selected molecules have similar Cα RMSD values with the positive control molecule (Supplementary Information, Figure S1). Ligand RMSDs were also investigated with two different cases: (i) Lig-Fit-Prot; (ii) LigFit-Lig. While Lig-Fit-Prot measures “translational” motion of the ligand at the binding pocket of protein, Lig-Fit-Lig measures “rotational” motion of the ligand. Heavy atoms of the ligands were considered in RMSD calculations. (Supplementary Information, Figure S2) shows Lig-Fit-Prot RMSDs. Vidarabine, xenazoic acid, tannic acid and mupirocin have RMSD values of 3.11, 4.69, 3.79 and 3.34 Å. Thus, they slightly change their initial docking poses. However, other two compounds (phytonodiol sodium diphosphate and cefpiramide) as well as reference molecule were structurally stable at the binding pocket. Corresponding values were found as 2.56, 2.58 and 1.98 Å for phytonodiol sodium diphosphate,

cefpiramide and

positive control molecule,

respectively. In the Lig-Fit-Lig ligand RMSDs all the selected ligands have RMSD values of less than 2.50 Å except the mupirocin, which may represent higher torsional motion of the molecule in the binding pocket of GzmB. The average Lig-Fit-Lig RMSD value was found as 2.03 Å for mupirocin. The corresponding values were 0.66, 1.52, 0.87, 1.45 and 1.71 Å for xenazoic acid, phytonodiol sodium diphosphate, vidarabine, tannic acid and cefpiramide, respectively. Reference molecule showed also small rotational motion throughout the simulations. Its average RMSD value throughout the simulation was found as 1.02 Å. (Supplementary Information, Figure S3) The RMSF is useful for characterizing local changes along the protein chain. The larger fluctuations were observed mainly at the loop regions, as expected (Supplementary Information, Figure S4). While reference molecule-bound complex structure has larger fluctuations (>5.0 Å) around the residue numbers of 56-68, in this region the maximum average RMSF values were around 2.0 Å for the studied compounds. Within the catalytic triad (His57, Asp102 and Ser195) region, none of the studied compounds showed larger fluctuations. Xenazoic acid and cefpiramide-bound complex structures showed slightly larger fluctuations (around 2.0 Å) near Asp102 residue (Supplementary Information, Figure S4).

16

Figure S5 at the Supplementary Information shows MM/GBSA results along 100-ns MD simulations. For this aim, 200-frames throughout 100-ns simulations were considered. The positive control molecule showed average MM/GBSA score of -86.00 kcal/mol. Within the selected hits, vidarabine showed not strong interactions compared to reference molecule and other selected hits. Cefpiramide showed the strongest interactions throughout the simulations with the target. Its average MM/GBSA score was found as 90.56 kcal/mol. The second best average score was found for the mupirocin which was -80.73 kcal/mol. The corresponding values for phytonadiol sodium diphosphate, tannic acid and xenazoic acid were -74.04, -68.03 and -77.06 kcal/mol, respectively. (Supplementary Information, Table S1) Since the studied structures have different molecular sizes (i.e., molecular weight), ligand efficiency scores (i.e., average MM/GBSA scores per heavy atoms) were also considered. While the ligand efficiency score of the reference molecule was found as -1.97 kcal/mol, corresponding values were -2.75, -2.30 and -2.15 kcal/mol for xenazoic acid, mupirocin and cefpiramide, respectively.

Figure 6. Changes of MD poses of Cefpiramide in the binding pocket of Granzyme B. Red, initial conformation; blue, final conformation.

Conformational changes of cefpiramide throughout the simulation are shown in a timescale coloring in Figure 6. In the figure, while initial binding pose was represented with

17

red color, blue color shows binding conformations of ligand at the end of the simulation. The more mobile region of the molecule throughout the simulations were found as pyridinone region at the binding pocket. Top-ranked conformation based on maximum occupancy within the reported binding pocket through hydrogen bonding and other potential non-covalent interactions were analyzed. To study the inhibition mechanisms of the ligands inside the binding pocket of GzmB, the top-6 molecules on the basis of MM/GBSA score and interaction with crucial residues were selected to investigate their interaction established with the active site amino acid residues. The selected molecules used in this study were; Tannic acid, Xenazoic acid, Cefpiramide, Vidarabine, Mupirocin, Phytonodiol sodium diphosphate and Reference molecule. In addition to the active site catalytic residues (His57, Asp102 and Ser195), other residues that are important in terms of their presence in the active site pocket are Lys40, Arg41, Gly216, Leu171, Arg226, Ser216, Gly193, Phe191, Lys192. Asp102 is not directly involved in catalytic reaction however it is used for charge stabilization during reaction mechanism. Cefpiramide made significant interactions with Ser100, Tyr174, Lys40, Arg226, Ser195, Lys192, His57, Gly216, Tyr215, that were present for 43%, 45%, 37%, 99%, 96%, 90%, 79%, 53%, 99%, of the simulation time via hydrogen bonds, hydrophobic interactions and π-cation interactions. Arg226 and Gly217 and Ser195 maintain 99% interactions throughout simulation time (Figure 7). Residues showing interactions in docking were also maintained during MD simulations. It is possible to have interaction >100% as some residues may have multiple interactions of a single type with the same atom. For example, the Arg side chain has three H-bond donors that can form hydrogen bonding interactions to a single H-bond acceptor as shown in Figure 7.

18

Figure 7: Bar graph showing interaction fractions (top-left) and interactions map during MD simulations for Cefpiramide (top-right); timeline representation (bottom) shows the interaction and contacts (H-bonds, Hydrophobic, Ionic, and water bridges). The top panel shows the total number of specific contacts the protein makes with the ligand over the course of the trajectory. The bottom panel shows which residues interact with the ligand in each trajectory frame. Some residues make more than one specific contacts with the ligand, which is represented by a darker shade of orange, according to the scale to the right of the plot.

For Mupirocin substantial interactions were retained with Ser214, Gly216 and Arg41 via hydrogen bonding interactions and water bridges (Supplementary Information, Figure S6).

Strong interactions were observed with Phytonodiol sodium diphosphate in docking. Most of the interactions were maintained during simulation time. Significant interactions

19

were made with His57, Lys40, Lys97, Lys192, Gly193, Asp194, Ser195, Arg226 that were existing for 66%, 40%, 99%, 66%, 31%, 50%, 88%, 100%, 99%, 91%, 76%, 82% of the simulation time through hydrogen bonds, some water bridges and π-cation interactions. Lys97 and Arg226 maintained 100 % interactions and made multiple hydrogen bonds (Figure 8).

Figure 8. 2D ligand interactions diagram for Phytonadiol sodium diphosphate complex throughout MD simulations. Tannic acid made ten hydrogen bonds in docking and other hydrophobic and non-polar interactions but only four hydrogen bonds were stable throughout MD simulations. A πcation interaction between the ligand and the Lys192 were stable in 41% of the simulation time. Constructed hydrogen bonds with Asn218 and Gly216 were stable and observed throughout the 37% and 42% of the simulation time, respectively. Corresponding values for the other hydrogen bonds were 44% and 37% for Tyr215 and Asp170, respectively. (Supplementary Information, Figure S7). Constructed chemical

20

interactions for the other selected compounds (Xenazoic acid, Vidarabine) were represented in Figures S8 and S9 at the Supplementary Information. Significant interactions were noted for the reference molecule. Tyr94, His57, Lys41 and Arg226 showed interaction throughout the simulation time (Supplementary Information, Figure S10). We also performed the replicas of each simulation for the selected molecules. MM/GBSA binding free energy values are re-calculated for the replica simulations and results can be found in Supplementary Table S1. As it can be seen in the table, all the molecules have similar binding free energies with first simulations results. For example, xenazoic acid has -77.06±7.34 kcal/mol predicted binding free energy values in the first MD simulations and in the replica simulations we calculated it as -74.32±7.37 kcal/mol. Indeed, all repeated simulations gave quite similar binding free energy values with the first run. Thus, it shows that obtained results are reliable. The RMSD and MM/GBSA plots of replicas can be found in supplementary materials Figures S11 and S12. All the ligands showed stable interaction with the active side residues as shown in simulation studies. These ligands can be used as a potential inhibitor for GzmB after required preclinical and clinical tests. As these are all FDA approved drugs, all these compounds are promising drugable compounds with therapeutic potential. All selected molecules had a higher docking score than reference molecule (-7.35 kcal/mol), proposing that selected molecules have the potential to inhibit the biological activity of GzmB. Docking score were higher for Tannic acid (-17.54 kcal/mol) followed by Mupirocin (-13.73 kcal/mol), Phytonodiol sodium diphosphate (-12.45 kcal/mol), Cefpiramide (-11.14 kcal/mol), Xenazoic acid (-10.75 kcal/mol), Vidarabine (-10.25 kcal/mol) and reference molecule (-7.35 kcal/mol). Glide SP and IFD algorithms were used for docking of all selected molecules and reference molecule as shown in Table 1. The top-docking results have been reported making hydrogen bond, hydrophobic interactions, polar and non-polar interactions with the catalytic residues and pocket specificity conferring residues of the GzmB (Figure 4). In the top-docking IFD pose, Cefpiramide showed hydrogen bonding with Tyr174, Gly193, Asn218, salt bridge with Arg226, some polar and hydrophobic interactions with Phe191, Lys215, Leu171, Tyr174,

21

175, and Phe99. Followed by Tannic acid showed similar interactions with Gly193, Arg41, Asp37, His57, Tyr174, π-π stacking interactions with Lys40 and some hydrophobic interaction with Tyr215, Leu171, Tyr174, Tyr175, Ala227, Phe99, Ile35, Leu39.

Figure 9: 2D and 3D ligand interactions diagram of top-docking pose of all selected molecules with GzmB. Hydrogen bonding, polar, nonpolar and hydrophobic interactions were noticed for topdocking pose of Mupirocin. In the corresponding docking poses, π-π stacking, π-cation, hydrogen bonding and salt bridge interactions are observed for Phytonodiol sodium

22

phosphate and Xenazoic acid. Vidarabine molecule showed hydrogen bonds with Arg41, Ser195 and Asn218, however salt bridge was formed with Arg226. Lys40 showed πcation interactions. (Figure 9). Reference molecule did not show interactions with catalytic triad in docking pose, however interactions involving hydrogen bonding, πcation, non-polar and hydrophobic interactions are shown for specificity conferring residues (Figure 9). For all the aforementioned complexes of GzmB, we have found that Tannic acid, Mupirocin, Cefpiramide, Xenazoic acid, Vidarabine and Phytonadiol sodium diphosphate can be used as a potential inhibitor of GzmB as they exhibit high predicted binding affinity (Table 1). It is concluded that the selected ligands were able to accommodate well inside the binding domain of GzmB by hydrogen bonds, salt bridges, hydrophobic, polar and nonpolar amino acids residues. We compared the scaffolds of 6 proposed molecules with reference molecule and with each other. Proposed 6 compounds have diverse scaffolds in their structures. However, this is not surprising since here we screen the large dataset of compounds and based on chemical interactions scores, a set of molecules were studied in detailed. Although the compounds may carry different structural fragments in their moieties, their nonbonded chemical interactions may have strong interactions, thus low binding energy (∆G) values. However, two of the selected compounds (cefpiramide and xenazoic acid) include peptide backbones similarly with the reference compound. Moreover, cefpiramide includes tetrazole moiety similarly with reference compound which has triazole fragment in its structure. When we check the structural similarities within the selected compounds, it is found that mupirocin and pytonadiol sodium diphposphate have long alkyl chain tails, vidarabine and pytonadiol sodium diphposphate compounds have phosphate fragments in their structures. MetaCore/ MetaDrug Analysis

GzmB is involved in various diseases like autoimmune diseases, cancer, heart failure, arthritis and chronic inflammation. We therefore, analyzed all the selected compounds for their potential therapeutic values and toxicities. This allowed us to get a thorough

23

understanding and prediction of their biological activity when compared to hundreds of other drugs using MetaCore/MetaDrug from Clarivate Analytics. [56, 57] With a cutoff value for therapeutic activity of 0.5 (predicted activity values are between 0 and 1) the selected hit molecules showed a significant therapeutic potential for heart failure, viral infections, arthritis, obesity, cancer, allergy, Schizophrenia and migraine. Our selected ligands have a higher prediction for their potential use for arthritis, cancer and heart failure (Table 3). Table

3: Therapeutic

activities selected

of FDA

approved molecules.

24

Toxicity test of reported GZMB inhibitors

Reference molecule and screened molecules were subjected to MetaCore/MetaDrug analysis for toxicity predictions [54]. In the toxicity profile check, we used 26-different toxicity models. Our selected molecules i.e. Cefpiramide and Xenazoic acids did not show any high toxicity predictions (Tables S2 and S3 in Supplementary Material). However, it is predicted that phytonadiol sodium diphosphate and tannic acid have high toxicity profiles in certain toxicity models such as cardiotoxicity. But selected molecules are FDA approved drugs, and they have passed through rigorous safety profiles, conferring their non-toxicity. Conclusions Limited data has shown the presence of inhibitors against granzyme B. For instance, PI-9 is the endogenous serine protease inhibitor (Serpin) that inhibits the activity of granzyme B inside cellular environment. Similarly, peptide substrate was also used for determining the substrate specificity of granzyme B however no detailed investigation about potential small molecules that has the ability of inhibit the activity of extracellular granzyme B is available. For a drug repurposing studies, all the inhibitors have to pass through rigorous tests/analysis i.e. therapeutic potential and toxicity tests that can then make them as potential candidates to be used as feasible inhibitors in clinical trials. Current drug repurposing study is a step ahead in finding a potential lead compound by screening database of FDA approved molecules. We have identified inhibitors (Tannic acid, Mupirocin, Phytonadiol sodium diphosphate, Cefpiramide, Xenazoic acid) with a potential of obstructing the active site of Granzyme B. Although, all identified inhibitors have shown stable interaction with the active site residues, based on the physiochemical

25

properties, Cefpiramide and Mupirocin can be considered as lead compounds for further drug repurposing studies.

References 1.

2.

3.

4. 5. 6.

7.

8.

9.

10. 11.

Ahmad, J., P.I. Bird, and D. Kaiserman, Analysis of the evolution of granule associated serine proteases of immune defence (GASPIDs) suggests a revised nomenclature. Biol Chem, 2014. 395(10): p. 1253-62. Caughey, G.H., A Pulmonary Perspective on GASPIDs: Granule-Associated Serine Peptidases of Immune Defense. Curr Respir Med Rev, 2006. 2(39): p. 263277. Baker, E., G.R. Sutherland, and M.J. Smyth, The gene encoding a human natural killer cell granule serine protease, Met-ase 1, maps to chromosome 19p13.3. Immunogenetics, 1994. 39(4): p. 294-5. Kaiserman, D., et al., The major human and mouse granzymes are structurally and functionally divergent. J Cell Biol, 2006. 175(4): p. 619-30. Peitsch, M.C. and J. Tschopp, Granzyme B. Methods Enzymol, 1994. 244: p. 807. Edwards, K.M., et al., The human cytotoxic T cell granule serine protease granzyme H has chymotrypsin-like (chymase) activity and is taken up into cytoplasmic vesicles reminiscent of granzyme B-containing endosomes. J Biol Chem, 1999. 274(43): p. 30468-73. Pilat, D., et al., The human Met-ase gene (GZMM): structure, sequence, and close physical linkage to the serine protease gene cluster on 19p13.3. Genomics, 1994. 24(3): p. 445-50. Plasman, K., et al., Conservation of the extended substrate specificity profiles among homologous granzymes across species. Mol Cell Proteomics, 2013. 12(10): p. 2921-34. Fu, Z., et al., Asp-ase Activity of the Opossum Granzyme B Supports the Role of Granzyme B as Part of Anti-Viral Immunity Already during Early Mammalian Evolution. PLoS One, 2016. 11(5): p. e0154886. Goping, I.S., et al., Granzyme B-induced apoptosis requires both direct caspase activation and relief of caspase inhibition. Immunity, 2003. 18(3): p. 355-65. Sutton, V.R., et al., Caspase activation by granzyme B is indirect, and caspase autoprocessing requires the release of proapoptotic mitochondrial factors. Immunity, 2003. 18(3): p. 319-29.

26

12. 13.

14. 15. 16. 17. 18. 19. 20.

21.

22.

23.

24. 25.

26.

27. 28.

Waterhouse, N.J., et al., A central role for Bid in granzyme B-induced apoptosis. J Biol Chem, 2005. 280(6): p. 4476-82. Andrade, F., et al., Granzyme H destroys the function of critical adenoviral proteins required for viral DNA replication and granzyme B inhibition. EMBO J, 2007. 26(8): p. 2148-57. Pardo, J., M.M. Simon, and C.J. Froelich, Granzyme A is a proinflammatory protease. Blood, 2009. 114(18): p. 3968; author reply 3969-70. van Eck, J.A., et al., A novel proinflammatory role for granzyme A. Cell Death Dis, 2017. 8(2): p. e2630. Metkar, S.S., et al., Human and mouse granzyme A induce a proinflammatory cytokine response. Immunity, 2008. 29(5): p. 720-33. Buzza, M.S. and P.I. Bird, Extracellular granzymes: current perspectives. Biol Chem, 2006. 387(7): p. 827-37. Joeckel, L.T., et al., Mouse granzyme K has pro-inflammatory potential. Cell Death Differ, 2011. 18(7): p. 1112-9. Kaiserman, D. and P.I. Bird, Control of granzymes by serpins. Cell Death Differ, 2010. 17(4): p. 586-95. Bird, C.H., et al., Selective regulation of apoptosis: the cytotoxic lymphocyte serpin proteinase inhibitor 9 protects against granzyme B-mediated apoptosis without perturbing the Fas cell death pathway. Mol Cell Biol, 1998. 18(11): p. 6387-98. Rotonda, J., et al., The three-dimensional structure of human granzyme B compared to caspase-3, key mediators of cell death with cleavage specificity for aspartic acid in P1. Chem Biol, 2001. 8(4): p. 357-68. Isaaz, S., et al., Serial killing by cytotoxic T lymphocytes: T cell receptor triggers degranulation, re‐filling of the lytic granules and secretion of lytic proteins via a non‐granule pathway. European journal of immunology, 1995. 25(4): p. 10711079. Jeng, M.Y., et al., Metabolic reprogramming of human CD8⁺ memory T cells through loss of SIRT1. The Journal of Experimental Medicine, 2018. 215(1): p. 51-62. Namekawa, T., et al., Functional subsets of CD4 T cells in rheumatoid synovitis. Arthritis & Rheumatism, 1998. 41(12): p. 2108-2116. Pardo, J., et al., Granzyme B is expressed in mouse mast cells in vivo and in vitro and causes delayed cell death independent of perforin. Cell Death And Differentiation, 2007. 14: p. 1768. Choy, J.C., et al., Granzyme B in atherosclerosis and transplant vascular disease: association with cell death and atherosclerotic disease severity. Modern pathology, 2003. 16(5): p. 460. Wagner, C., et al., Granzyme B and perforin: constitutive expression in human polymorphonuclear neutrophils. Blood, 2004. 103(3): p. 1099-1104. Tschopp, C.M., et al., Granzyme B, a novel mediator of allergic inflammation: its induction and release in blood basophils and human asthma. Blood, 2006. 108(7): p. 2290-2299.

27

29.

30. 31.

32. 33. 34.

35. 36.

37. 38. 39.

40.

41.

42.

43.

44. 45.

Rissoan, M.-C., et al., Subtractive hybridization reveals the expression of immunoglobulinlike transcript 7, Eph-B1, granzyme B, and 3 novel transcripts in human plasmacytoid dendritic cells. Blood, 2002. 100(9): p. 3295-3303. Turner, C.T., D. Lim, and D.J. Granville, Granzyme B in skin inflammation and disease. Matrix Biology, 2017. Korkmaz, B., et al., Neutrophil elastase, proteinase 3, and cathepsin G as therapeutic targets in human diseases. Pharmacological reviews, 2010. 62(4): p. 726-759. Hendel, A., et al., Granzymes in age-related cardiovascular and pulmonary diseases. Cell death and differentiation, 2010. 17(4): p. 596. Granville, D., Granzymes in disease: bench to bedside. 2010, Nature Publishing Group. Shen, Y., et al., Topical small molecule granzyme B inhibitor improves remodeling in a murine model of impaired burn wound healing. Experimental & molecular medicine, 2018. 50(5): p. 68. Casciola-Rosen, L., et al., Granzyme B: evidence for a role in the origin of myasthenia gravis. Journal of neuroimmunology, 2008. 201: p. 33-40. Afonina, I.S., et al., Granzyme B-dependent proteolysis acts as a switch to enhance the proinflammatory activity of IL-1α. Molecular cell, 2011. 44(2): p. 265-278. Omoto, Y., et al., Granzyme B is a novel interleukin-18 converting enzyme. Journal of dermatological science, 2010. 59(2): p. 129-135. Wang, T., et al., Granzyme B-induced neurotoxicity is mediated via activation of PAR-1 receptor and Kv1. 3 channel. PloS one, 2012. 7(8): p. e43950. Ikemoto, T., et al., Plasma granzyme B as a predicting factor of coronary artery disease—clinical significance in patients with chronic renal failure. Journal of cardiology, 2009. 54(3): p. 409-415. Buzza, M.S., et al., Extracellular matrix remodeling by human granzyme B via cleavage of vitronectin, fibronectin, and laminin. Journal of Biological Chemistry, 2005. 280(25): p. 23549-23558. Lee, R.T., et al., Circumferential stress and matrix metalloproteinase 1 in human coronary atherosclerosis: implications for plaque rupture. Arteriosclerosis, thrombosis, and vascular biology, 1996. 16(8): p. 1070-1073. Chamberlain, C.M. and D.J. Granville, The role of Granzyme B in atheromatous diseases. Canadian Journal of Physiology and Pharmacology, 2007. 85(1): p. 8995. Sastry, G.M., et al., Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments. J Comput Aided Mol Des, 2013. 27(3): p. 221-34. Bas, D.C., D.M. Rogers, and J.H. Jensen, Very fast prediction and rationalization of pKa values for protein-ligand complexes. Proteins, 2008. 73(3): p. 765-83. Jorgensen WL, Tirado-Rives J (1988). The OPLS Force Field for Proteins. Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 110 (6): 1657–1666

28

46.

Shelly JC, Cholleti A, Frye LL, Greenwood JR, Timlin MR, Uchimaya M. Epik: a software program for pK(a) prediction and protonation state generation for drug-like molecules. J. Comput Aided Mol Des 2007, 21(12): 681-91

47.

Estebanez-Perpina, E., et al., Crystal structure of the caspase activator human granzyme B, a proteinase highly specific for an Asp-P1 residue. Biol Chem, 2000. 381(12): p. 1203-14. Friesner, R.A., et al., Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. J Med Chem, 2004. 47(7): p. 1739-49. Friesner, R.A., et al., Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J Med Chem, 2006. 49(21): p. 6177-96. Desmond Molecular Dynamics System, D. E. Shaw Research, New York, NY. Hoover, W.G., Canonical dynamics: Equilibrium phase-space distributions. Phys Rev A Gen Phys, 1985. 31(3): p. 1695-1697. Martyna, G.J., Remarks on "Constant-temperature molecular dynamics with momentum conservation". Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics, 1994. 50(4): p. 3234-3236. Li, J., et al., The VSGB 2.0 model: a next generation energy model for high resolution protein structure modeling. Proteins, 2011. 79(10): p. 2794-812. MetaCore/MetaDrug, https://portal.genego.com

48.

49.

50. 51. 52.

53. 54. 55. 56.

57.

Kim, M.S., et al., Approaches to design non-covalent inhibitors for human granzyme B (hGrB). Org Biomol Chem, 2014. 12(44): p. 8952-65. Matthews, E.J., et al., Assessment of the health effects of chemicals in humans: I. QSAR estimation of the maximum recommended therapeutic dose (MRTD) and no effect level (NOEL) of organic chemicals based on clinical trial data. Curr Drug Discov Technol, 2004. 1(1): p. 61-76. Matthews, E.J., et al., Assessment of the health effects of chemicals in humans: II. Construction of an adverse effects database for QSAR modeling. Curr Drug Discov Technol, 2004. 1(4): p. 243-54.

29

Highlights •

Granzyme B is a cytotoxic serine protease having the potential to initiate apoptosis in a cancerous and virally infected cells.

•

FDA approved drugs were screened in virtual screening to find potential inhibitors against Granzyme B.

•

Docking of top compounds showed strong interactions with the active site catalytic triad and other specificity conferring residues.

•

Top six compounds were further used for long MD simulations.

We have no conflict of interest to declare. Prof Dr Serdar DURDAGI