Structure Optimization of Neuraminidase Inhibitors as Potential Anti-Influenza (H1N1Inhibitors) Agents Using QSAR and Molecular Docking Studies

Introduction

Influenza virus infection is commonly known as flu and it is the contagious etiologic agent that causes an acute respiratory infection; hence it has always been a major threat (1) to human health worldwide and cause for economic burden. Currently, effective chemotherapy for influenza virus is also limited due to newly discovered drug resistance in mutant strain (2). Although the viral replicative cycle has revealed several potential molecular targets (M2 proteins endonuclease, hemagglutin and neuraminidase) that can be used as anti-flu for drug design (3-6). There are currently only a few licensed drugs available for influenza treatment. M2 inhibitors such as Amantadine and Rimantadine, which act specifically against Influenza A virus by blocking the ion channel of the M2 protein, provide only limited protection due to a narrow spectrum of activity (7). Neuraminidase (NA), also called sialidase, is the major surface glycoprotein that possesses enzymatic activity essential for viral replication and infection (8). The crystal structure of influenza virus Neuraminidase is well known and it is a critical enzyme in the infestation, replication, maturity and delivery of influenza virus (9).

Though Neuraminidase inhibitors like Zanamivir displays excellent anti-viral activity when administered intranasal, it is less effective when delivered systemically. It has very low oral bioavailability and is rapidly eliminated by renal excretion (10).Oseltamivir is orally active, but it has been reported to cause vomiting, nausea and several allergic reactions. While remaining Neuraminidase inhibitors such as Peramivir, Laninamivir are still in phase III clinical trials. Peramivir shows less oral bioavailability than Oseltamivir. So, there is still an enormous need to design and identify new agents for the chemotherapy of influenza virus infection and formulate effective drugs for systemic administration. Therefore in present study we put forth following objectives:

1. To search the structural fragments responsible for selective NA inhibition.

2. To quantify the contribution of each structural fragments towards biological activity by comparing.

3. The contributing physicochemical parameters of required chemical groups in the given series.

The understanding of the enzyme’s structure and functions can be helpful for the development of novel Neuraminidase Inhibitors. Quantitative Structure Activity Relationships (QSAR) for different sets of compounds has been reported by Verma and Hansch (11) which gives seventeen different QSAR equations to understand chemical and biological interactions governing their activities toward influenza neuraminidase. Another QSAR model was developed with spatial, topological, electronic, thermodynamic and E-state indices on 30 thiourea analogues as influenza neuraminidase inhibitors using SYBYL 7.1 molecular modelling software and genetic algorithm method (12).

In this study, we performed 3D-QSAR using kNN-MFA method with the help of V-Life MDS 3.5. Considering electronic, steric parameters stated by 3D-QSAR studies, and also considering structural necessities for selective neuraminidase inhibitors at active binding sites of neuraminidase, the pharmacophore of acyl thiourea was optimized so as to inhibit neuraminidase efficiently. Optimized pharmacophore was used and little bit modified to generate basic templates so as to design new chemical entities using Combi-Lib tool of V-Life MDS 3.5. Generated library contains more than 100 entities, designed on the basis of equation of multiple linear regression studies. Out of 100, best 10 molecules were selected on the basis of predicted activity. Those best 10 molecules were subjected to docking studies. Docking studies were carried out by Glide tool of Schrodinger Inc. Software. QSAR (13-17) model gave good insights for developing new analogues as influenza virus NAI. On the basis of docking studies, it is concluded that designed entities show best binding interactions including Hydrogen bonds as well as Van-der-Waals forces within neuraminidase active binding site so as to inhibit neuraminidase efficiently for the treatment of influenza.

Experimental

Dataset

A series (18) of total 28 compounds for whose absolute IC₅₀ values is reported for their anti-influenza activity was used for QSAR studies (Table 1 and Table 2). Biological activity was expressed in terms of pIC₅₀ which was calculated as log1/IC₅₀ against neuraminidase enzyme. All QSAR models were generated using a Training set of 22 molecules. Test set of 6 molecules with varied chemical and biological activities were used to access the predictive power of generated QSAR models using training set of molecules. The sphere exclusion method was used for the selection of molecules in training and test sets.

Uni-Column statistics for training set and test set were generated to check correctness of selection criteria for trainings and test set molecules (Table 1).Selection of molecules in the training set and test is a key and important feature of any QSAR model, therefore due care was taken in such a way that biological activities of all compounds in test set lie within the maximum and minimum value range of biological activities of training set of compounds.

The common pharmacophore used for QSAR study is Acyl thiourea with R₁ and R₂ substitutions is represented (Figure 1).

Acyl thiourea template

Figure 1

Selected series of acyl thiourea compounds with biological activity as neuraminidase inhibitors is depicted in tabular form along with training set and test set molecules (Table 2). Selected series is a collection of structures having wide range of neuraminidase inhibitory activity ranging from minimum 0.08 to maximum 18.5 µM.Validation of QSAR studies: Models generated by 3D-QSAR study was cross-validated using standard LOO procedure (19).

QSAR Studies

3D-QSAR Studies

3D-QSAR studies were performed using V-Life Molecular Design Suite Software Version 3.5 (20). 3D-QSAR Studies were carried out by k Nearest Neighbor Molecular Field Analysis (kNN MFA) using Simulated Annealing (SA) variable selection method majority of its k-Nearest Neighbors in the training set (21-22).

3D model was generated using following steps:

1. Molecules were optimized MMFF energy minimization method before alignment. Optimization is necessary process for proper alignment of molecules around template.

2. kNN-MFA method requires suitable alignment of given set of molecules; alignment was carried out by template based alignment method. Alignment of Substituted Acyl Thiourea Derivatives using template based alignment method is represented below (Figure 2).

1. This was followed by generation of common rectangular grid around the molecules, the steric and electrostatic interaction energies were computed at the lattice points of the grid using a methyl probe of charge +1.

2. The optimal training and test set were generated.

3. Models were generated by SA kNN MFA method.

4. The activity of the test set of compounds was predicted.

5. The 3D-QSAR model was evaluated using important statistical measures such as, ‘n’ is number of molecules (> 20 molecules), ‘k’ is number of descriptors in a model (statistically n/5 descriptors in a model), ‘q^2’’ is cross-validated r² (>.0.5), ‘pred_r²’ is r² for external test set (>.0.5), ‘q²_Se’ is Standard Error in cross validation and ‘pred_r²_se’ is Standard Error in external validation.

Alignment of Substituted Acyl Thiourea Derivatives using template based alignment method.

Figure 2

Optimization of Pharmacophore

The information obtained from 3D QSAR study and active binding site was used to optimize the electrostatic and steric requirements as well as essential structural requirements for neuraminidase inhibition around the Acyl Thiourea nucleus for selective inhibition of Neuraminidase and in turn to enhance anti-influenza activity.

Generated 3D descriptors are correlated structurally and biologically with most potent compound (Figure 3) from the selected series of acyl thiourea having neuraminidase inhibitory activity up to 0.08 µM.

Most potent Compound from selected series

Figure 3

Considering above requirements obtained from QSAR studies and binding site studies, pharmacophore was optimized.

Design of new chemical entities using designing of new templates and generation of library

Information obtanined from QSAR studies and active binding site have helped us a lot to design new chemical Entities (NCE‘s).We have generated two templates using Lead Grow tool of V-Life MDS software and on the basis of these templates,we have generated more than hundrded molecules using CombiLib tool of V-Life MDS software.

Molecular Docking studies

Molecular docking studies were performed using Glide (5.0) (26) tool of Schrodinger molecular docking software. Those who had predicted activity as similar as pIC50 of most active molecule present in the reported series were selected for docking studies.

ADMET Prediction

All designed compounds were filtered by predicting their Absorption, Distribution, Metabolism, Excretion and Toxicity (ADMET) properties by means of Qikprop 2.2 Tool of, Schrodinger (Table 8).

Results and Discussion

Result of 3D-QSAR

Series of Acyl thiourea derivatives was carefully divided into training set and test set. 3D-QSAR Studies were carried out by k Nearest Neighbor Molecular Field Analysis (kNN MFA) on the selective training set and test set as depicted in Table 2.

kNN MFA is a novel method which utilizes k-Nearest Neighbor (kNN) principle to correlate molecular field descriptors with biological activity. The kNN methodology relies upon a simple distance learning approach. In this method an unknown member is classified according to the majority of its k-Nearest Neighbors in the training set. The nearness is measured by an appropriate distance metrics (e.g., a molecular similarity measure calculated using field interactions of molecular structures). The standard kNN Method is implemented simply as follows:

1. The distances between an unknown object (u) and all other objects in the training set were calculated.

2. The k objects were selected from the training set most similar to object u, according to the calculated distances.

3. The object u was classified with the group to which the majority of the k objects belong. An optimal k-value is selected by optimization through the classification of a test set of samples or by Leave-One Out (LOO) cross validation (19). The variables and optimal k-values were chosen using different variable selection methods. Here we have used simulated annealing as variable selection method.

Statistical result of 3D-QSAR kNN MFA method is tabulated below (Table 3).

Stereo view of the molecular rectangular field grid generated around the superposed molecular units of Acyl Thiourea series using SA kNN-MFA Model

Figure 4

Model showed best internal as well as external predictivity as q² = 0.6840, pred_r² = 0.7697 and also error occurred during internal and external validation obtained was very low as q²_se = 0.2811, pred_r²_se = 0.2972.Stereo view of superposed molecular units of Acyl Thiourea series shows the generated data points around the pharmacophore in three dimensional manners, which is shown below (Figure 5).

Stereo view of the molecular rectangular field grid generated around the superposed molecular units of Acyl Thiourea series using SA kNN-MFA Model

Figure 5

Thus, SA KNN-MFA model leads to identification of various local interacting molecular features responsible for activity variation and hence provides direction for design of new molecules in a convenient way.

Interpretation from 3D QSAR Studies and Correlation of 3D descriptors

Electronic Parameters

3D QSAR studies revealed the electronic requirements around the acyl thiourea pharmacophore. The points those were found in SA KNN-MFA model are E_449, E_450 implying that these points are indeed significant for structure activity relationship and require the electronic properties as mentioned in the ranges in parenthesis for maximum biological activity.

Range for electronic descriptor E_450 (-0.0456, 0.0118).Descriptor ranging from negative to positive values specifies that the pharmacophore should be substituted with electronegative and electropositive groups so that it should increase the neuraminidase inhibiting activity, such as –OH, –Cl, -NO₂, -OR and most potent compound possesses two electronegative groups such as –Cl, -OC₂H₅ from the parent series.

Steric Parameters

The steric points that were found in SA KNN-MFA model is S_333 implying that these points are indeed significant for structure activity relationship and require the steric properties as mentioned in the ranges in parenthesis for maximum biological activity. Range for steric descriptor is S_333 (-0.0053, -0.0041). Descriptor showing negative values specifies that the pharmacophore should be substituted with sterically less bulky group so that it can increase the NA inhibiting activity, such as lower alkyl groups, unsubstituted or mono substituted aryl rings.

Study of Structural Necessity for Neuraminidase Inhibition

Wang and Wade have reported guidelines for structural modification of developing neuraminidase inhibitors and the design of novel inhibitors in order to optimize inhibitory activity (Figure 6) (23). Since the catalytic site of influenza virus Neuraminidase was totally conserved among all influenza viral strains, an ideal drug which effectively blocked one neuraminidase would be effective at blocking all other neuraminidases, even those on viruses which have not yet appeared in humans (8).

Thus, along with QSAR study, we are herein considering the essential structural necessities which should be present in the acyl thiourea pharmacophore. It will be beneficial in such a way that we can design better compounds which would fit into neuraminidase binding site perfectly and would show best interactions than standard, marketed neuraminidase inhibitors.

Important structural features for potent inhibitor and corresponding neuraminidase residues (23)

Figure 6

A negatively charged (Nc) group that makes strong charge-charge interactions with the triarginyl pocket as Arg 292, Arg 371 and Arg 118 is highly favorable for binding. Positively charged (Pc) group is favoured for binding by the electrostatic interactions with residue Asp151.While, Arg 152 requires hydrogen bond acceptor (Ha) group for binding. Trp 178 requires such a structural functionality which should be hydrogen bond donor (Hd) as well as small hydrophobic group for favourable binding. While, large hydrophobic group contributing to the van der Waals interaction should be preferred for three different amino acid residues such as Arg 224, Ala 246 and Glu 276.

Along with structural necessities for neuraminidase inhibition, we are considering air plane model of neuraminidase enzyme, so as to design potential derivatives with excellent anti-influenza activity (24). Wang and co-workers (24) derived an airplane model of the neuraminidase active site as illustrated in (Figure 7) to summarize the basic structural requirements of potent neuraminidase inhibitors. The active site of neuraminidase has four main well conserved binding sites. The centre of site 2 is about 6 A.U. from site 1 and about 4 A.U. from site 3, while site 4 is about 6 A.U. from site 1 and 5 A.U. from site 3. Sites 1 and 3 are separated by 9-10 A.U. or about 7 single bond lengths.

Diagram of neuraminidase sites S1-S4 and important nearby residues (25).

Figure 7

Gong and co-workers have proposed about the pocket study of neuraminidase active site (8).

By taken into consideration, structural necessity and airplane model and pockets of neuraminidase enzyme, we will design such new chemical entities which are much potential than standard and selective inhibitors.

Design of New Chemical Entities

Designed compounds generated in this way were then screened by three type of screening methods; Lipinski’s rule and prediction of activity using multiple linear regression equation [pIC₅₀ = 0.250765 T_C_N_6 + 0.106694 T_C_O_2 + 0.637398 MMFF_2 - 0.623115 chiV3 + 6.1896] obtained by the two dimensional QSAR studies and to ensure drug like pharmacokinetic profile by prediction of ADMET properties for finding a new compounds having neuraminidase inhibitor and anti-influenza activity (20). New chemical entities were generated using optimized pharmacophore shown in (Figure 8).

Pharmacophore Requirement of Acyl thiourea Derivatives Generated for Selective Inhibition of neuraminidase for anti-influenza activity

Figure 8

More than one hundreds of molecules were generated using CombiLib tool which follows the Lipinski’s rule, but we have selected only 10 most active molecules on the basis of predicted activity. The most potent compounds have positive values of pIC₅₀, similarly the predicted activity values for new chemical entities were found to lye towards positive side.Compounds qualifying all required parameters set for Lipinski’s Screen/filter are indicated by ADRXWS strings. The columns containing the Lipinski’s screen score and other column containing the strings of alphabets, ADRXWS indicate that all 6 conditions are satisfied by that corresponding compound. Lesser the screen score, lesser is the pharmacokinetic compatibility (drug likeliness) for that designed compound.

Generation of Templates

We have generated two templates considering above stated pharmacophore requirements. Designed templates consist of dicarbonyl groups which are different than parent series, while parent series has single acyl group in it, we have placed another acyl group instead of any aryl ring containing oxygen atom, as in parent series has furan ring in it (19).

One template (Figure 9) consists of 3 –NH groups in the pharmacophore, while another template (Figure 10) consists of 2 –NH groups, one –NH group is isosterically replaced by –CH₂ group.

Template A

Figure 9

Template B

Figure 10

Both of the templates have seven single bonds between R¹ and R² substituents as it is basic structural requirement for neuraminidase inhibition, where in case of R¹ and R², we have used aryl ring with electronegative and electropositive groups such as–OH,-COOH,-NO₂ and –OCH_3, –NH₂ and aryl ring also contributes to hydrophobicity, as correlated in 3D-QSAR studies. Those seven single bonds are represented in red colour in designed template (Figure 11). Even placement of extra carbonyl functionality and aryl ring with electronegative carboxyl group in the acyl thiourea pharmacophore becomes helpful for binding it with triarginyl residues as discussed above (Figure 6,7).

Designed template with active binding site

Figure 11

Generation of Combinatorial Library

We have generated more than hundrded molecules using CombiLib tool (20) which follows the Lipinski’s rule, but we have selected only 10 most active molecules on the basis of their predicted activity using Multiple Linear regression equation.

We have selected only four molecules from the B template (Table 4) considering promising predicted activity.While from A template we have selected seven molecules (Table 5) on the basis of predicted activity.

Out of total 11 designed new chemical entities only 3 compounds viz. R-8,R-9 and R-10 show predicted activity less potent as compared with most potent compound from the parent series with pIC50 1.09691,because higher the value of pIC50 more potent is the entity.

Molecular Docking studies

Molecular Docking is a key tool in structural molecular biology and computer-assisted drug design. Selected most active molecules were docked on crystallographic structure of neuraminidase enzyme available in the RCSB PDB Database (Code: 3b7e) co-crystallized with the ligand Zanamivir (27). A molecular docking study helps to determine possible interaction of new chemical entities with the enzyme on PDB (3b7e).

The goal of ligand–protein docking is to predict the predominant binding mode(s) of a ligand with a protein of known three-dimensional structure. Virtual screening on the basis of molecular descriptors and physicochemical properties of active ligands has great usefulness in finding hits and leads through library enrichment for screening (28), a strategy that is also well-used for reducing and enriching the library of ligands for molecular docking; there are recent reports that ligand shape-matching does as well as, if not better than, docking (29).

Glide was found to produce least number of inaccurate poses and 85% of Glides binding models had an RMSD of 1.4 A⁰ or less from native co-crystallized structures (30). The Glide docking program approximated a complete systematic search of the conformational, orientation and positional space of the docked ligand molecules in to the receptor (protein) binding pocket.

Overview of docking of New Chemical Entities

All the designed compounds that show good predicted activity and follows Lipinski’s rule were docked into neuraminidase enzyme (pdb code: 3b7e) to study the binding mode of designed compounds. Further screening to sort out the best compound having good binding affinity which was compared with binding mode of Standard Neuraminidase Inhibitors, e.g. Oseltamivir, Zanamivir and Peramivir results of which are depicted in (Table 6 and 7).

The reliability of the docking results was first checked by comparing the best docking poses obtained for the co-crystallized inhibitor with its bound conformation.

As a result, a root mean square deviation (RMSD) of 0.7 Å was found suggesting that the docking procedure could be relied on to predict the binding mode of our compounds.

Poses of interactions of ligands along with standard are represented in (Figure 12, 13, 14 and 15).

Hydrogen bond interaction of R-7

Figure 12

Hydrogen bond interaction of R-10

Figure 13

Hydrogen bond interaction of Oseltamivir

Figure 14

Hydrogen bond interaction of Peramivir

Figure 15

Key findings of overall Docking studies

The close inspection of result of molecular docking studies indicated that the designed compounds docked better than ligand Oseltamivir and Peramivir. Following are the reasons for the same.

1. G-Score

Molecules R-7 and R-8 show highest G-score than both the standards Peramivir and oseltamivir. While, molecules such as R-2, R-4, R-5, R-9 and R-10 show higher G-score than oseltamivir.

2. Hydrogen bond interactions

Molecule R-10 shows highest hydrogen bond interactions i.e. 8 than Peramivir and Oseltamivir, which show 7 and 4 interactions respectively. While, R-6 shows 6 hydrogen bond interactions which are more than oseltamivir.

3. Van der Waals interactions

R-1, R-3, R-4, R-5, R-6 show less bad van der Waals interactions than oseltamivir and Peramivir, where only Peramivir shows highest bad van der Waals interactions as 14 and oseltamivir shows 8 bad van der Waals interactions. Molecules R-1, R-3 and R-6 show zero ugly forces, where oseltamivir and Peramivir show 2 and 1, respectively.

4. In all the bad and ugly contacts penalize the G-score and overall conformation of ligand does matter a lot while binding to amino acid residues at active binding site. This also penalizes the energy of the model.

Key interactions are depicted in tabular form (Table 7). It was observed from docking studies that all ligands lie in the same pocket of neuraminidase binding site of enzyme containing Arg 152, Tyr 406, Arg 371, Arg 118, Glu 227, Glu 277, Asp 151, Arg 292, Asn 294 and Ser 246 amino acids. Thus designed compounds showed a good binding affinity for interaction with neuraminidase active binding site.R-7, R-8, R-9 and R-10 show Hydrogen bond interaction with different arginine residues, in this case one carbonyl group, a negatively charged group from dicarbonyl thiourea interacts with positively charged groups or ions such as –NH₂ or H⁺ of OH arginine residues. Here, one of structural necessity for selective neuraminidase inhibition is satisfied.

R-7, R-8, R-9 and R-10 show electrostatic interaction with different glutamic acid residues. In their case, either a positively charged group such as –NH₂ or positively charged ion H⁺ of –OH or –COOH interacts with carbonyl or carboxyl group of glutamic acid residues, satisfying structural necessity for selective neuraminidase inhibitor.

Oseltamivir shows Hydrogen bonding with amino acids Arg 152, Tyr 406, Arg 371, Arg 118 which are one of essential interactions for neuraminidase inhibition, which is described previously in study of structural necessity for neuraminidase inhibition. Distance of hydrogen bonds in NCEs and standard compounds also explain better interaction (Table 7) with neuraminidase binding site in 3b7e pdb.

The hydrogen bond distance in case of compound R-7 with Glu 227 (1.793 A.U.) which is less as compared to the hydrogen bonding interaction in case of Peramivir which is 2.477 A.U. Hence, less hydrogen bonding distance than standards is required for better binding which is reported in NCEs during docking study.

G-score is calculated taking all above aspects in to consideration, while designing new chemical entities to inhibit Neuraminidase more effectively and in turn to optimize the pharmacophore required for selective inhibition of Neuraminidase.

ADMET Prediction

Prediction of ADMET properties was used as last screen to sort out those compounds that already follow Lipinski’s rule and show good predicted activity and binding conformation at Neuraminidase receptor.

The parameters illustrated in (Table 8) Qikprop analysis show significant results. CNS parameter is related with absorption of entity through Blood brain barrier, standard limit for CNS is -2 to +2, where -2 shows inactive CNS penetration and +2 shows active CNS penetration. All the designed entities show satisfactory results, with negative values, indicating poor CNS penetration. % Oral Absorption parameter is related with extent of oral absorption of drug, indicating suitable route of administration, if it is going to be formulated. If entity shows more than 80% oral absorption, it is considered to be highly absorbed. While if any entity shows less than 25% oral absorption, it is considered to be poorly absorbed. Metabolites suggest the number of metabolites which will possibly generate after undergoing metabolic changes, number of metabolites should range from 1-8.

Conclusion

The thorough analysis of results of 3D QSAR studies and the structural necessities for neuraminidase inhibition have helped us to make a decision about the electronic, steric, hydrophobic nature of substitution pattern around the selected acyl thiourea pharmacophore. Two different Templates were generated on the basis of above information using lead grow tool of V-life Molecular design suite. Using the information, the New Chemical Entities (NCEs) were designed using CombiLib tool with the help of multiple linear regression equation and activities were also computed for the designed NCEs.

The major reason for failure of NCEs at latter stages of drug discovery process i.e. drug like pharmacokinetic profile set up using Qikprop (33) 2.2 Tool of, Schrodinger; so that only drug like NCEs would be generated and resultant NCEs would not have the pharmacokinetic inadequacies. The generated NCEs were analyzed by Lipinski’s screen. Results indicated that designed NCEs are satisfying all the parameters set for Lipinski’s screen. The most potent derivatives were subjected to molecular docking studies to get further insights of interactions of NCEs with neuraminidase.

In present study, we performed molecular modeling study demonstrating that acyl thiourea derivatives inhibit neuraminidase by binding to the site. The results of dry lab work will be analyzed thoroughly to find out correctness of the rational used for the design of NCEs in general and optimization of pharmacophore for selective neuraminidase inhibition of in particular.