Design of new potent HTLV-1 protease inhibitors: in silico study

Introduction

Human T lymphotropic virus type 1 (HTLV-1), a member of Delta-retrovirus family, has single-strand RNA genome and positive polarity [1]. This virus is known to cause a T-cell lymphoma cancer directly and is endemic to various locations around the world, including Japan and Caribbean Sea. It is estimated that approximately 15-20 million people worldwide are infected with HTLV-1 [2-4]. HTLV-1 infects human blood cells and causes a variety of diseases in humans, such as adult T-Cell Leukemia (ATL), HTLV-1 associated myelopathy /tropical spastic paraparesis (HAM/TSP) and neurological disorders [2-4]. HTLV-1 protease, an aspartic protease essential for maturation and replication cycle of HTLV-1 virus, is a homodimeric protein with each chain consisting of 125 residues and is highly specific to its substrate [5, 6]. Two acidic residues are located in the narrow tunnel-shaped active site of enzyme accommodating substrates or inhibitors [7]. This protein is responsible for the processing of Gag and Gag-pro-pol polyprotein during virus maturation and catalyzes an essential step in virus replication cycle [6].

Many compounds have been developed to specifically inhibit protease activity such as statine-based inhibitors; however, they have modest potential against proteases. Most of these compounds are only active at micromolar concentrations [8]. Furthermore, small organic molecules such as drugs have high bioavailability and are resistant to proteolitic degradation [9]. Nowadays, peptidomimetics, compounds that act as appropriate substitutes for receptor interacting peptides, have attracted many pharmaceutical scientists [10, 11]. These molecules are typically designed based on modification of an existing peptide, or similar systems that mimic peptides. These compounds are more beneficial due to properties such as high stability and biological activity [12, 13]. New series of inhibitors can be developed by modifying the structures of another class of peptidomimetics that are peptoids and peptidosulfonamides [14]. Therefore, design of the K… N… I… (KNI) series of ligands based on peptidomimetics is appealing strategy since these ligands display strong inhibitory effects [8]. HIV protease is an aspartic protease as HTLV protease [6, 15]. The sequence similarity of these two proteins is 28% while the similarity in the active site is 45%. [16, 17].

In this work, novel protease inhibitors were designed against HTLV-1 protease. Most of them were designed with a sulfonamide or peptoid functional group in their backbone. First, absorption, distribution, metabolism and excretion (ADME) properties of ligands were estimated according to Lipinski's rule and the toxicity and carcinogenicity of new inhibitors were evaluated. Subsequently, these ligands were docked to HTLV-1 protease. The docking calculations offered best ligands to continue. Some of the previous synthesized inhibitors (KNI series) also were docked to HTLV-1 protease in order to compare their inhibition constants with our designed ligands. Finally, SP-4 and SP-5 that showed best properties as protease inhibitors were selected to estimate complex stability. It should be noted that the specific inhibitors for HTLV-1 protease also have the ability to inhibit HIV protease but the specific inhibitors for HIV protease cannot inhibit HTLV-1 protease [8, 18]. To investigate inhibitory effect, SP-4 and SP-5 were also docked to HIV-protease. Our results showed that the novel designed ligands (SP-4 and SP-5) strongly inhibit HTLV-1 and HIV-proteases. Our expectation is that the obtained results open new horizons to in human T-cell lymphoma virus treatment.

Materials and methods

Ligand and receptor preparation: Ligands were constructed using HyperChem program [20] (Fig 1). Then a solvation box created for each ligand. Molecular dynamics simulation was performed at 300 K with 0.5 picoseconds runtime. Periodic boundary conditions were performed with step size 0.005 ps until molecules obtained their most energetically favored structures [20].

ADME properties of the compounds were predicted using molsoft web services (Molsoft LLC., unpublished, http://molsoft.com/mprop) which estimates the absorption, distribution, metabolism and excretion of candidate compounds. The percentage of Drug-Likeness for all designed drugs was also calculated by this server. The Lipinski's rule is considered as a simple way to estimate drug-likeness [21]. Toxicity and carcinogenicity of new inhibitors, as well as their effects on organisms, were evaluated using the web service LAZAR Toxicity Predictions [22].

The structure of HTLV-1 protease with inhibition KNI-10681 inhibitor (PDB ID|3LIT) and HIV protease (PDB ID| 4HE9) were obtained from PDB database [7, 23]. Inhibitors and all solvent molecules were removed from the protein structure.

Molecular Docking: In this study, Autodock was used as a docking software; the grid point spacing was set to 0.375 Å and cubic box with 70 × 70 × 70 dimensions. All docking parameters were calculated using Lamarckian Genetic Algorithm (LGA). All selected amide bonds were non-rotatable and docking was iterated one hundred [24, 25]. Two top-scoring compounds were selected, both to confirm the docking result and find the best complex.

Molecular dynamics simulations: Molecular dynamics simulations were also performed on complexes that showed highest inhibitor constant level at docking. Moreover, dynamic behaviors of inhibitors in the active site were assessed. The simulation was performed using Gromacs version 4.6.3 package and energy was calculated by Gromos96 force field parameter set 53a6 [26]. The docking output from HTLV-1 protease (PDB ID | 3LIT) was used by Gromacs. All-atom MD simulations were also performed. Two sets of simulations, with similar parameters and setup, were performed for HTLV-1 protease in the presence of SP-4 and SP-5 compounds, separately. HTLV-1 protease was solvated with explicit solvent (SPC water) [27] in a cubic box, which left 1.0 nm space around the solute for HTLV-1 protease. Subsequently, the net charge of system was neutralized. After energy minimization with the steepest descent method, the system was equilibrated by a 100 ps simulation. It should be noted that during this time interval, positions of all bonds were constrained. The system was coupled to a Berenson’s thermostat method with the reference temperature at 300 K [28]. Neighbor searching and long-range electrostatics was done using PME [29, 30]. This was followed by full scale MD simulation. Coordinates were saved every 2 ps to the trajectory. MD simulation was performed during 20ns for all minimized structures for each of four simulations with 2 fs time steps.

Results AND Discussion

To determine the compounds from designed inhibitors, we scored them according to LAZER, Molsoft (unpublished, link) and Lipinski’s rules server. The compounds are shown in Table 1 and 2. The SM nomenclature is arbitrary/abbreviation for Sulfona Peptoid. From this point on, we refer to the determined compounds as SM-1, SM-2, SM-3, SP-4 and SP-5. LAZAR program suggested that there was no evidence of carcinogenic activity in mouse and rat while all ligands had carcinogenic effect on hamster. The molecular docking is frequently used in pharmaceutical research to predict and promote appropriate compounds for rational drug design [26]. Theoretically, all fifteen molecules showed very good binding energy ranging from -8.89 to -12.16 kcal/mol (Table 3). The inhibition constant values for SM-1, SM-2, SM-3 and SP-4 compounds were 268.43, 305.58, 202.03 and 47.69 nM, respectively. To compare the accuracy of our result, KNI series ligands were docked into the HTLV-1 protease using same condition (Table 3). Based on results, the new designed inhibitors significantly improve inhibition constant compared to KNI series ligands in HTLV-1 complex. The binding energy for compounds SP-4 and SP-5 was -9.99 kcal/mol and -12.16 kcal/mol with inhibition constant 47.69 and 1.23 nM, respectively. The logarithm of ligand partition coefficient (logP) is 3.65 for SP-4 and 2.5 for SP-5. Their Drug-Likeness score ranges from 0.09 to 1.69. Among them the Drug-Likeness scores for SP-4 and SP-5 is 0.87 and 1.07, with approximate molecular weight of 572.27 and 558.29, aqueous solubility (logS) value of -6.40 and -4.94, Polar surface area (PSA) value of 117.20 and 126.37, nine similar hydrogen bond acceptors and finally, five and nine hydrogen bond donors, respectively. It should be noted that the maximum Drug-Likeness calculated by the Molsoft software is 1. PSA describes drug permeability calculated by Molsoft. The drugs with a PSA greater than 140 Å are unable to enter into the cell membrane and a PSA less than 60 Å is necessary to penetrate the blood-brain barrier. Thus SP-4 and SP-5 are potentially interesting for further optimization with molecular dynamic simulation. It is postulated that eight subsites belonging to the two symmetrical chains of the HTLV-1 protease form that active site. These subsites are annotated as S1, S2, S3, S4, and S1’, S2’, S3’, S4’ depending on the chain they belong to. The residues located in the HTLV-1 active site include Arg10, Lys95, Asn96, and Asn97in S3 subsite;Asp36, Met37, Asn53, Thr54, Ser55, Cys90, and Val92 in S4 subsite; Leu30, Gly34, Val56, Leu57, Gln96, Gln97, Trp98 in S1subsite, and two catalytic aspartyl residues (Asp32) positioned in both symmetrical chains of protease [16, 31]. According to the docking results, the strongest H-bond interaction has occurred between SP-4 and SP-5 inhibitors with both S3 and S4 subsite in first chain of HTLV-1 protease.

RMSD results during Molecular Dynamic simulation illustrated that the absolute structure of SP-5 and SP-4 and their complex with HTLV-1 protease were stable during running time (Fig 2).

Conversely, patterns of hydrogen bonds on the ligand-protease complex fluctuate during MD simulation. An extra H-bond was formed between Arg10 in S3 subsite in the first chain of HTLV-1 protease and the semi-circular form of SP-4 whereas H-bond between SP-5-Gly 34 was weakened and a new H-bond was formed with Ala 59 after MD simulation running time interval (Fig 3). It seems that the conformational rearrangements of ligands after MD running time may cause the alteration of hydrogen bond number between ligands and HTLV-1 protease. Electrostatic interactions in the SP-5-HTLV-1 protease complex were also increased.

According to the literature, most of synthesized or designed inhibitors for HTLV-1 protease are located as a linear form which is typical for peptidic and peptidomimetic ligands bound to retroviral proteases, on a cavity formed between two symmetrical chains of protease through H-bonds and electrostatic interactions [8, 32]. Therefore, insertion of cyclic-aromatic hydrocarbons to both N-terminal and C-terminal of new designed ligands affected hydrophobic interactions between ligand and protease. These hydrophobic forces caused the ligand conformation conversion from linear to semi-circular form. It is worth mentioning that KNI-1595b, one of the synthetic inhibitors introduced by Satoh et al., containing two aromatic rings in both ends of ligand indicates better inhibition constant (21.9 nM) in KNI inhibitor series [8].

The subsites S3/S3' in HTLV-1 protease have three unique residues (Leu57', Asn97, Trp98) and three residues identical to other retroviral enzymes (Arg10, Leu30, Asp36) [29]. Based on the data in the literature for KNI inhibitor series, protease structure remains unchanged in the presence of different ligands whereas the loop containing residues 96-98 (loop 96-98) exhibit relative variability due to different functional group on ligands [8].  Loop96-98 plays a significant role in protease activity, although it is located outside the active site [8]. Li Mi et al. suggested that the presence of the unique Trp-98 and Leu-57' significantly changes the nature of S3/S3' pockets in HTLV-1 protease and make it hydrophobe. They further showed that Asn97 in loop 96-98 is shielded by the large side chain of Trp98 and does not interact with linear form inhibitors [32, 33]. The two aromatic rings in both C-terminus and N-terminus of SP-4 or SP-5 semi-circular inhibitors surround the large side chain of Trp98. The t-phenyl group in C-terminal of KNI-series is linked to protease through hydrophobic interactions with hydrophobic side chain of alanine, valine, leucine and phenylalanine residues. These hydrophobic interactions prevent any conformational change at the C-terminus of inhibitor [8]. Hydrophobic interactions occur between the aromatic planes in phenyl group at SP-4 or SP-5 N-terminus, Trp98 and naphthalene group at SP-4 or SP-5C-terminus (Fig 4). It seems that these hydrophobic interactions reduceTrp98 flexibility as well as loop96-98 and desheild Asn-97 to interact with inhibitor and especially have accommodated Trp98 in second chain of HTLV-1 protease.HIV protease is a homodimer protein similarly as HTLV protease, with each subunit composed of 99 amino acids [33,34]. The active site lies between the identical subunits and has the characteristic Asp-Thr-Gly (Asp25, Thr26 and Gly27) sequence common to aspartic proteases. The two Asp25 residues (one from each chain) act as the catalytic residues. The active site of HIV protease consists of Leu23,Pro81,Gly27' residues in S1^’ subsite; Asp30' Asp29', Ala28', Gly48', Ile50 in S2^’ subsite; Asp25', Ile84', Gly49, Pro81', Gly27, Arg6' in S1 subsite and Asp30, Asp29, Val32, Gly48, Ile50' in S2 subsite.

To investigate whether new designed inhibitor can inhibit HIV-protease, SP-4 and SP-5 were also docked to HIV-protease. Based on docking results, SP-4 and SP-5 compounds inhibit HIV-protease with binding energy -9.99 kcal/mol and -8.78 kcal/mol and inhibition constant 47.53 and 376.56 nM, respectively. Molecular dynamic simulations were performed for both of them at 20ns. RMSD results during Molecular Dynamic simulation illustrated that the absolute structure SP-5 and SP-4 HIV-protease complexes were stable during running time (Fig 6).

Patterns of hydrogen bonds on ligand- protease complex fluctuate during MD simulation. It seems that the conformational rearrangements of ligand, especially for SP-4 after, MD running time may cause alteration of hydrogen bond number between ligands and HIV protease. The pattern of hydrogen bonds was changed in SP-4-HIV protease complex. In this case, however, the residues involved in H-bond shifted from Arg8 and Asp25 to Ile50 and Asp29 in S1 subsite. Moreover, SP-4 interacted with Asp 25 residue of HIV protease electrostatically (Fig7). The SP-5 compound also connected to Asp25 and Asp25’ and strongly inhibit HIV-protease.

It is concluded that hydrophobic force plays an important role in suppressing protease activity. It is also proposed that designing and development of novel compounds based on aromatic hydrocarbons in both terminus of inhibitors and sulfanamid groups very promising for efficient treatment. These results suggest that insertion of naphthalene to the N-terminus of Sulfonapeptoid inhibitors as a new functional group makes them suitable ligands to inhibit HTLV-1 protease. Furthermore, we propose that the designed inhibitors (SP-4 and SP-5) with the proper modification have the ability to be developed as a novel pharmaceutical compound.

Conflict of Interest: The authors declare no conflict of interest.