Molecular insights of protein contour recognition with ligand pharmacophoric sites through combinatorial library design and MD simulation in validating HTLV-1 PR inhibitors

Abstract

Retroviruses HIV-1 and HTLV-1 are chiefly considered to be the most dangerous pathogens in Homo sapiens. These two viruses have structurally unique protease (PR) enzymes, which are having common function of its replication mechanism. Though HIV PR drugs failed to inhibit HTLV-1 infections, they emphatically emphasise the need for designing new lead compounds against HTLV-1 PR. Therefore, we tried to understand the binding level interactions through the charge environment present in both ligand and protein active sites. The domino effect illustrates that libraries of purvalanol-A are attuned to fill allosteric binding site of HTLV-1 PR through molecular recognition and shows proper binding of ligand pharmacophoric features in receptor contours. Our screening evaluates seven compounds from purvalanol-A libraries, and these compounds’ pharmacophore searches for an appropriate place in the binding site and it places well according to respective receptor contour surfaces. Thus our result provides a platform for the progress of more effective compounds, which are better in free energy calculation, molecular docking, ADME and molecular dynamics studies. Finally, this research provided novel chemical scaffolds for HTLV-1 drug discovery.

---

Introduction

The retrovirus human T-cell lymphotropic virus type 1 (HTLV-1) was discovered in the early 1980s. HTLV-1 infection is associated with the development of adult T-cell leukemia/lymphoma (ATL/ATLL), a clonal aggressive malignancy of CD4⁺ T lymphocytes.^1,2 At present, HTLV-1 infects approximately 20 million individuals all around the world and it is the first retrovirus to be unambiguously linked causally to a human cancer.³ In resemblance to human immunodeficiency virus (HIV), HTLV-1 mainly infects CD4 T-cells, which are the central regulators of the acquired immune response.⁴ To establish a persistent infection, HTLV-1 perturbs the regulation of CD4 T cells, which sometimes leads to adult T-cell leukemia (ATL).^5,6 When compared to other retroviruses, HTLV-1 encodes the protease (PR) enzyme, which is essential for viral maturation.⁷ Therefore targeting the protease enzyme HTLV-1 PR prevents viral proliferation and maturation, which makes the PR enzyme a key drug target for the development of new potential lead compounds.⁸ Tozser et al. (2007) stated that the mechanism by which the virus causes diseases is still not known but their studies indicated that viral replication is critical for the development of HTLV-1 associated myelopathy, and initial studies suggested that blocking replication with PR inhibitors had a therapeutic effect.⁹ Therefore, based on the success of HIV-1 PR inhibitors, HTLV-1 PR is also considered as a potential target for chemotherapy.¹⁰ But recent research reports state that anti-HIV PR drugs cannot function as HTLV-1 PR blockers and several successful HIV-1 PR inhibitors failed to provide the inhibitory activity against HTLV-1 PR,^8,11 so prediction of a potent inhibitor for HTLV-1 PR is highly essential for human welfare.

Li et al. (2005) crystallographically solved the X-ray structure of a HTLV-1 PR and reported that the HTLV-1 PR enzyme is an attractive drug target for anti-cancer drug design.¹² Structure of HTLV-1 PR represents the homodimer, containing 125 residues per chain, which is a longer sequence compared with HIV-1 PR.¹³

The two chains of the HTLV-1 PR monomer are bound by non-bonding interactions with the active site at the interface between two monomers.¹⁴ Similar to HIV PR, the HTLV-1 PR also shows the same position of the active site between the two monomeric chains, but HTLV-1 PR is elongated by an extra 10 amino acid residue “tail” at the C-terminus.¹⁵ Therefore the active site region is expanded in HTLV-1 PR, and the known HIV PR drugs cannot sustain a inhibitory profile against HTLV-1 PR.¹⁶ The structural information of HTLV-1 PR complexed with an inhibitor (PDB id 3LIN, 3LIQ, 3LIT, 3LIV, 3LIX and 3LIY) shows that active sites are located between the two monomer chains comprising of ARG10, LEU30, ASP32, GLY34, ALA35, ASP36, MET37, VAL39, LEU56, LEU57, ALA59, LEU91, TRP98 and LLE100 amino acids. These residues are known to contribute to the binding interactions and have the capacity to accommodate a ligand and inhibit the drug target.¹²

The presence of active site information boosts structure-based drug discovery to design HTLV-1 PR inhibitors, while improved understanding of active site contour features is an important reference for pocket specific de novo chemical library screening and scoring functions.¹⁷ A combined approach of receptor environment and screening of libraries with protein–ligand charge distribution was computed by the scoring function to recognize the importance of pharmacophoric features, which elucidates the arrangement of chemical features that are shared by molecules exhibiting activity at a receptor.^18–20 Currently, predicting accurate binding free energies of new leads in particular receptor sites is challenging and from this work, we have consistently introduced a number of structural adaptation strategies on the original structure of purvalanol-A, which induces viral proliferation and apoptosis.^21,22 Here, charge calculation of a receptor active site for predicting chemical requirement (contour maps) was carried out using SiteMap, which has the potential to accept the ligand with respective pharmacophoric features (chemical sharing of receptor–ligand). Thus the contour maps based interaction analysis between ligand and protein is crucial for developing new inhibitors. The protein contour maps are charged environments that interact only with their respective pharmacophoric site of small molecules.^23–25 The hydrogen bond acceptor, hydrogen bond donor, hydrophobic group, negatively charged group, positively charged group and aromatic ring pharmacophoric features in ligands tend to bind with structural factors of HTLV-1 PR surface of active site macromolecule as shown in Fig. 1 and Fig. S1a and b (ESI†). Contour-based lead optimization enhances the vitality of understanding to find out effective and potent drug for HTLV-1 PR.^25–27

Fig. 1 The contour based ligand design from the structure of HTLV-1 protease, chain monomer subunits are shown as red and grey ribbons and contours are analyzed for ligand rearrangement with pharmacophore rearrangement.

Materials and methods

System configuration

All studies were carried out in a high-performance cluster operated with Cent OS V6.2 Linux operating platform. The hardware specifications are HPC cluster-Super micro SC826TQ-R1200 LIB series, running with two Intel Xeon E5620 Quad Core 2.4 GHz processors on 32 GB RAM. The software specifications included are the academic version of Gromacs v4.5 for (molecular dynamics) MD simulations and the commercial version of Schrodinger 2012 software package (Version 9.2; LLC, New York, NY) for docking protocol.

Molecular modeling environmental setup

The typical structure file available from the PDB was not suitable for immediate use in molecular modeling calculations and so the crystal structure of HTLV-1 PR (PDB id: 2B7F) was prepared through protein preparation wizard.¹² The two monomers chain A and B were taken for subsequent investigation. By using the Prime, missing residues and the missing loops were filled from the SEQRES records in the PDB file. Firstly, the bond orders were assigned, hydrogen atoms were added and all the crystallographic water molecules without any contacts were removed. Minimization was performed until the average RMSD of the non-hydrogen atoms to reach 0.3 Å. The default sampling and water optimization using OPLS-2005 force field and impref minimization were applied to prepare the apo protein. The reported inhibitor of HTLV-1 PR, namely purvalanol-A, was prepared with Lig-Prep 2.4 using OPLS-2005 force field. The preparation was done so as to retain the original state of ligand and chirality. Up to 32 conformations per ligand structure were generated using the default energy ring confirmation.²⁸

Binding site and contour map prediction

Druggability sites were located through clustering the favorable regions by means of vdW charges on protein surface using SiteMap 2.4.²³ It implements OPLS-2005 force field parameters to estimate the interaction energies of probes placed at all points on three dimensional grids that encompass the entire protein. To identify the top ranked potential druggability sites, it required at least 10 site points, and then the environment was set to be more restrictive by definitions of hydrophobicity, using standardized grid and crop site maps at 3 Å.²⁹ The requirement of physico-chemical properties in ligands to design the best inhibitor against the protein active site contours were predicted by using SiteMap.³⁰ The docked results of known inhibitors were subjected to physico-chemical properties of protein contours prediction. The ligand molecule was then picked manually to evaluate the single binding region around the inhibitor and an additional region of around 6 Å buffer for examination.²³

Combinatorial library design

Based on protein physico-chemical requirements on the active site, the rearrangement regions were noted in the purvalanol-A structure and by using the ligand designer script R¹, R², R³ and R⁴ sites were picked (Fig. S1a and b, ESI†).³¹ Addition of atoms (chain structure, cycle structure, C-groups, miscellaneous, protecting groups details are given in ESI,† Fig. S2) from ChemSketch program³² and of radicals in each R-position was made and imported in Maestro and optimized through OPLS-2005 force field using ligprep.²⁸

Molecular docking simulation

The interactions between the protein and ligand were carried out by the molecular docking program Glide V.6.1.³³ Here two different type of docking methodology were performed. For purvalanol-A, experimental active site based docking was performed and for the designed inhibitors druggability site based docking was performed.^34–37

Receptor information based grid generation

For purvalanol-A, the receptor grid generation was processed manually by picking the ligand entry and specifying the centroid of specific residues including ARG10, LEU30, ASP32, GLY34, ALA35, ASP36, MET37, VAL39, LEU56, LEU57, ALA59, LEU91, TRP98, and LLE100.¹² The position of spherical region that should be occupied by a particular ligand during docking was set on XYZ axes with co-ordinates of 92.80, 53.33 and 54.16, respectively. The same active site residues were picked for the HBONDS constraints and were designated as flexible residues.³⁷

Site-based grid generation

The druggability site based on the top-ranked Site Score in SiteMap was prearranged as Glide input files and the receptor grid generation was carried out with the white colored spheres along the protein.³⁸ To recognize the ligand position, the entry of white spheres was picked using the van der Waals radius scaling factor of 1.0 and partial charge cutoff of 0.25. Here the position of the grid box was set on XYZ axes with the measurements of 87.94, 55.59, and 54.43, respectively, with radius 2.0.³⁹

Druggability site and receptor based ligand docking

Both receptor and druggability site based docking was performed with XP docking protocol.³³ Here, Glide generated conformations internally and passed those conformations through a series of filters. In the first stage, the ligand center at various grid positions of 1.00 Å was placed and the ligands were allowed to rotate around the three angles. While at this point of time dissimilar binding modes of the ligands were removed based on the crude score values and geometrical filter. In the second filter stage a grid-based force field evaluation and refinement of docking solutions including torsional and rigid body movements of the ligand was analyzed through the OPLS-AA force field. A small number of surviving docking solutions were subjected to a Monte Carlo procedure for minimized energy score. Subsequently the final energy evaluation was done with the Glide Score and a single best pose was generated as output for a particular ligand with the help of following equation:

G-score = a × vdW + b × Coul + Lipo + Hbond + Metal + BuryP + RotB + Site

where vdW = van der Waals energy; Coul = Coulomb energy; Lipo = lipophilic contact term; HBond = hydrogen-bonding term; Metal = metal-binding term; BuryP = penalty for buried polar group; RotB = penalty for freezing rotatable bonds; Site = polar interaction at active site; and the coefficient of vdW and Coul are: a = 0.065, b = 0.0130.

The above equation gives the G-score, which includes the main factors of van der Waals energy, Coulomb energy, hydrophobic grid potential, hydrogen-bonding term, metal-binding term, buried polar groups, freezing rotatable bonds and polar interactions in the active site.^34,35

Physiochemical descriptor calculation

The eMBrAcE, Prime MM-GBSA, Liaison and Qikprop calculations for the ligand–receptor complex structures were performed to generate the ligand and structure-based descriptors (LSBD).⁴⁰ eMBrAcE applies multiple minimizations using OPLS-2005 with the surface generalized born implicit water solvent model. Constant dielectric electrostatic treatment were applied specifically 1.0 (default) and non-bonded cutoff a criterion is extended and iterations are counted up to 5000, which calculates molecular mechanics energy minimization of the complex.^41,42 The permeability, adsorption, dissolution, metabolism, and excretion (ADME) properties of compounds were done with Qikprop.⁴³

Molecular recognition with pharmacophores

The best compounds from docking and free energy analysis were selected to create pharmacophore sites (site points) using PHASE. PHASE supplies a built-in set of six pharmacophore features, which includes hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic group (H), negatively charged group (N), positively charged group (P) and aromatic ring (R).^44,45 Pharmacophoric features are defined by a set of chemical structure patterns that form the active part of the drug molecule, particularly for structure based targets. Based on these pharmacophoric sites, the binding contours of the active site are visualized for molecular recognition.⁴⁶

Molecular dynamics simulation

The MD simulation studies were carried out for the crystal structure of HTLV-1 PR on a 20 ns timescale in order to understand the stability and intra-molecular conformational changes occurs in the protein. The GROMACS program package (http://www.gromacs.org) adopting the OPLS-AA force field parameters was used for energy minimization and MD simulations.⁴⁷ For the MD simulation studies, the structure was solvated using the TIP3P water model, and the solvated structure was energy minimized using steepest descent method, terminating when maximum force was found to be smaller than 100 kJ mol⁻¹ nm⁻¹.⁴⁸ The total simulation was performed in the NPT ensemble at constant temperature (300 K) and pressure (1 bar), with a time step of 2 fs. Meanwhile, NVT was performed for 1 ns, and the minimized structure was equilibrated with a timescale of 20 ns.^49,50 Additionally, MD simulations were performed for the ligand bound docked structures of HTLV-1 PR, which were the outcome of combinatorial library design and free energy analysis. The initial structure of the receptor and ligands were cleaned using GROMOS96 force field and then the topology files were generated for the receptor and ligands separately using PRODRG tool.⁵¹ The simulation system was created manually by importing the ligand topology into the system pursued along with a dodecahedron box with a margin of 1 nm and the system was filled with water using the SPC explicit solvation model. The system was applied with energy minimization and the atomic velocities were adjusted according to Maxwell Boltzmann distribution at 300 K with a periodic scaling of 0.1 ps. A presimulation run of 20 ps was applied to relax the system and to remove the geometric restrains which eventually appeared at the initialization of the run. All the simulations were carried out at constant pressure and temperature (NPT) ensemble. The Berendsen coupling was employed to maintain a constant temperature of 300 K and constant semi-isotropic pressure of 1 bar with coupling time of 2.0 fs and the coordinates were saved. The simulation timescale for the ligand bound-form was 20 ns and RMSD analysis was performed for understanding the stability of ligands.⁵²

Results and discussion

Active site of HTLV-1 PR – experimental and theoretical analysis

The comparison of experimental and theoretical active site will enhance the binding mode and also cause involvement of a few other amino acids as active sites. The information of the active site in both experimental (crystallographic information) and theoretical prediction provides strong support for the computational analysis of protein/small molecule interactions. In this study based on the available complex-bound HTLV-1 PR crystal structures with statine inhibitors, we understand that amino acids of ARG10, LEU30, ASP32, GLY34, ALA35, ASP36, MET37, VAL39, LEU56, LEU57, ALA59, LEU91, TRP98 and LLE100 were designated as the active sites for the crystal structure of human T-cell leukemia virus protease (PDB id: 2B7F). Theoretically predicted druggability regions coding amino acids are observed in maestro, and show white colored spheres appearing between two monomer chains in the Z-shaped surface (Fig. S3, ESI†). Through the knowledge based method, the druggability sites are noted as ARG10, LEU30, LEU31, ASP32, THR33, GLY34, ALA35, ASP36, MET37, THR38, VAL39, SER 55, LEU56, LEU57, GLY58, ALA59, LEU91, ASN97, TRP98, LLE100. But meanwhile we noticed that a few amino acids in the theoretical study are not experimentally proven as active sites, and so molecular docking of library compounds may suggest the presence of new active sites. In our previous works, we have mentioned that the residue Met37 has a mutational effect and so, here the Purvalanol-A has redocked with M37D residue (Fig. S4, ESI†). Interestingly, the wild type interaction residues do not show the interactions role with M37D residue and Purvalanol-A has lost its ability to bind with HTLV-1 PR. The reported compounds from this work, if they have the ability to bind with Met37, will also have the ability to bind with mutant HTLV-1 PR.^53,54

Purvalanol-A interaction and generation of contour maps

Binding of purvalanol-A to HIV-1 PR has not yet reported, but in case of HTLV-1 it has been reported that Purvalanol-A has the effect of inhibition.²² To understand interaction of the known HTLV-1 PR active purvalanol-A inhibitor, we have docked purvalanol-A with known active site residues. It shows very good interactions and bound complex between both chain A (GLY34) and chain B (ASP32, GLY34) amino acids (Fig. 2). These molecular interactions are considered for contour map generation to check physio-chemical property requirements of receptor active site. Through this contour map analysis, the whole active site area of the receptor protein is invigilated and their structural properties are elucidated. Depending on the results, we understand that HTLV-1 PR requires plenty of hydrophilic region (HPR) from the ligand molecules; hydrogen bond acceptor regions (HBAR), hydrogen bond donor regions (HBDR) and hydrophobic regions (HPBR) are available in range; and absence of metal binding regions (MBR) in receptor active site. The surface area of each region of the contour maps with the respect to purvalanol-A is reported in Table 1. In visualization of contour maps, each HPR (green mesh), HBAR (red mesh), HBDR (blue mesh) and HPBR (yellow mesh) are shown in Fig. 3.

Fig. 2 The reported compounds purvalanol-A showing interactions with both chain A (GLY34) and chain B (ASP32, GLY34) amino acids.

Table 1 Represents the molecular interaction between the crystal structure of HTLV-PR and known inhibitor purvalanol-A

Molecular docking in HTLV-PR (2B7F)		Contour maps prediction of purvalanol-A
Ligand name	Purvalanol-A	H-bond acceptor region (HBAR)	738.927
G-score	−6.702560	H-bond donor region (HBDR)	1102.740
No. of interactions	3(A1 and B2)	Hydrophilic region (HPR)	1803.663
Interacting residues	GLY34, ASP 32	Hydrophobic region (HPBR)	311.759
Chain A	GLY34	Metal binding region (MBR)	0.000
Chain B	GLY34, ASP 32

---

Fig. 3 Visualization of contour maps showing the active site environment of HTLV PR protein, which includes hydrogen bond acceptor region (HBAR – red mesh), hydrogen bond donor region (HBDR – blue mesh), hydrophilic region (HPR – green mesh) and hydrophobic region (HPBR – yellow mesh) in the presence of purvalanol-A.

Combinatorial library design

The ligand rearrangement position is marked with placement of receptor contours and R¹, R², R³ and R⁴ positions are drawn (Fig. S1, ESI†) and based on this, more replacements are done at the R positions (derivatives of purvalanol-A R¹, R², R³, R⁴ are provided in the ESI,† Fig. S2). The glide docking program was employed as a primary docking search engine to dock library compounds with that predicted druggability regions of HTLV-1 PR. The white colored sphere covers the 4 Å area of a binding pocket, representing a binding cleft where the ligand molecule can bind specifically at the region. When ligand molecules replaced the white color dots in the binding cleft region, the atomic structure arrangements attained their electronic level and sharing of hydrogen bonds between the acceptor and donor of the ligands and receptors occurred. Resulting docking studies were evaluated for the binding mechanism of library compounds with HTLV-1 PR and had a more negative value of G-score, which indicated good binding affinity of the ligand with the receptor (ESI,† Tables S1–S4). On assessment with all library compounds, the best compounds were chosen, which came top in the G-score criteria and form better interactions. In this obtained conformation of ligands, the numbers of hydrogen bonds formed are noted with respective amino acid involvement.

Rationalization of active sites in ligand recognition

Theoretical active site prediction using SiteMap predicts more correlation in results with experimental active sites; confirmation of active sites are theoretically checked by libraries of purvalanol-A interaction results with better G-scores than known purvalanol-A and are considered for active site residue analysis. In a refined docking environment, hydrogen bonding interactions are fully based on donor and acceptors present in the protein and ligands. The atoms in the donor and acceptor regions form the charge-based interactions, i.e., the positive and negative regions between the hydrogen bond donor (HBD) and hydrogen bond acceptor (HBA) are able to interact. The atoms with the capacity to donate their hydrogen atom towards the acceptor regions to form bonding bridges are called HBD atoms and the atoms with a tendency to accept the donated hydrogen atoms are called HBA atoms. The hydrogen bonding formed between the HBD and HBA throughout the ligand–receptor interaction was accounted and active sites were analyzed (Table 2). Each interaction is tabulated and provided in Tables S1–S4 of the ESI.† The libraries of purvalanol-A that show better results than the source inhibitor in terms of G-score are analyzed with active site information for site matching for ligand patterns with reference experimental active site information. Based on G-score and energy parameters, the nine best compounds are chosen, i.e. fluroenylmethoxycarbonyl (R¹ – A), phosphate (R² – B), phosphate (R⁴ – C), benzyloxycarbonyl (R⁴ – D), phosphate (R² – E), sulfate (R³ – F), phenyloxycarbonyl (R⁴ – G), trifluroacetyl (R³ – H) and fluroenylmethoxycarbonyl (R⁴ – I) (Fig. 4). These nine compounds show good interactions with HTLV-1 PR, with good scoring and are well bound inside the binding pocket (Table 3). The better results in the molecular docking inform about the interacting active site residues and the comparison of experimentally available and theoretically predicted results is cross-validated with docking results in specification of residue capacity to bind ligand atoms. Theoretical prediction of active sites may or may not be accurate, so comparison of theoretical prediction with reference to experimentally available and site based molecular interaction results with reporting residues having capacity to interact was performed. Comparison of active sites of experimentally reported theoretically predicted and docking interaction residues results with LEU56, LEU91 and LLE100 shows the absence of active site property and residues SER55, GLY58 and ASN97 are identified as having active site property, showing these residues can also function as active site of HTLV-1 PR.

Table 2 The active site difference in both experimental and theoretical prediction, which is from reported PDB complex 3LIN, 3LIQ, 3LIT, 3LIV, 3LIX, 3LIY, crystallographic information of active site is obtained and theoretical prediction is obtained from SiteMap

Experimental active sites reported		Theoretical active site prediction		Docking interaction residues		Additional residues interaction
ARG10	LEU30	ARG10	LEU30	ARG10	LEU30	SER55
ASP32	GLY34	LEU31	ASP32	ASP32	GLY34	GLY58
ALA35	ASP36	THR33	GLY34	ALA35	ASP36	ASN97
MET37	VAL39	ALA35	ASP36	MET37	VAL39
LEU56	LEU57	MET37	THR38	SER55	LEU57
ALA59	LEU91	VAL39	SER55	GLY58	ALA59
TRP98	LLE100	LEU56	LEU57	ASN97	TRP98
		GLY58	ALA59
		LEU91	ASN97
		TRP98	LLE100

---

Fig. 4 Molecular docking interactions of nine best lead compounds passed in both docking score and binding energy calculations (A = R¹ fluroenylmethoxycarbonyl, B = R² phosphito, C = R⁴ phosphato, D = R⁴ benzyloxycarbonyl, E = R² phosphato, F = R³ sulfato, G = R⁴ phenyloxycarbonyl, H = R³ trifluroacetyl, I = R⁴ fluroenylmethoxycarbonyl. These nine compounds are having good G-score, binding energy, activity, Hbond interaction.)

Table 3 The best ligand extracted from molecular docking, binding energy descriptors, and pharmacophore features of the ligand, which are denoted as A, D, H, N, P and R

Ligand name	Pharmacophore feature						Computed values of protein–ligand interaction
	A	D	H	N	P	R	Gscore	DGbind
Fluroenylmethoxycarbonyl (R^1 – A)	6	3	4	0	0	5	−8.565992	−42.570530
Phosphate (R^2 – B)	7	5	3	0	0	3	−8.590773	−48.911172
Phosphate (R^4 – C)	8	3	3	1	0	3	−8.520319	−40.469488
Benzyloxycarbonyl (R^4 – D)	6	3	3	0	0	4	−9.017948	−48.872633
Phosphate (R^2 – E)	8	3	3	1	0	3	−8.711677	−44.911172
Sulfate (R^3 – F)	6	2	4	1	0	3	−8.468488	−43.204161
Phenyloxycarbonyl (R^4 – G)	6	3	3	0	0	4	−8.812543	−39.448403
Trifluroacetyl (R^3 – H)	5	2	5	0	0	3	−9.141660	−49.290872
Fluroenylmethoxycarbonyl (R^4 – I)	6	3	4	0	0	5	−9.221705	−50.243652

---

Ligand and structure based descriptors

LSBDs are generated to maximize the capabilities of the methods for predicting and rank-ordering the binding affinities of compounds for a given target protein. Here, ligand–receptor complex descriptors are obtained from the structural information and it was calculated by Embrace, Liaison, Prime MM/GBSA, and Qikprop, which are tabulated and used for enhancing the docking results. Based on these results, binding free energy was calculated for each docking pose. From that, compounds which have tendency to execute more than −35 kcal mol⁻¹ are considered for next phase of study. The good G-score, binding energy, and H-bond interaction of these limited screened compounds are tabulated in Table 3. The nine screened compounds have a good enough Gscore, energy parameters and also have tendency to form strong interactions with the HTLV-1 PR. When comparing these nine compounds with available inhibitors, each compound has unique pharmacophoric sites, but all compounds are dominated by hydrogen bond acceptors, hydrophilic regions and aromatic rings.

For these nine compounds, various energy parameters has been evaluated, such as electrostatic energy (U_elec), van der Waals energy (U_vdW) and solvation energy (G_solv) using Embrace and Liaison. These are energetics of ligands binding to active site surfaces and are able to be used to determine the binding mode of known high-affinity ligands and find new active compounds. The solvent water is used for calculation of binding energy; specifically, U_cav is a prediction of cavity energy, by using the RADAP (apparent radius of the solvent) and vdW energy is predicted by using the RAD_HS (hard sphere radius of the solvent). These solvent-based energy approaches are physical terms, which are contributing in protein ligand molecular recognition that computationally predicts the energy parameters using hybrid Monte Carlo simulation with dielectric constant. The energy values in Table 4 show the energy configuration between the ligand and protein binding; these energies are responsible for binding of new inhibitors with HTLV-1 PR. ADME delivers a medicinal chemistry task – PSA, represented in Table 4. This is polar surface area, determing the induction of new lead capacity into the cell membrane; the new molecules of C and E having PSA pf more than 140 are eliminated by this medicinal chemistry task.

Table 4 Represents OPLS-2005 based vdW energy, electrostatic energy, solvation energy in ligand binding and cavity energy in presence of hybrid water model as solvent through hybrid Monte Carlo simulation approach and polar surface area of selected compounds

Computed values of energies involvement in present of hybrid water model
Ligand	vdW energy	Electrostatic energy	Solvation energy	Liaison (UvdW)	Liaison (Ucav)	Liaison (Uelec)	PSA
A	−261.210000	−100.380000	185.720000	−81.129000	9.360600	3.667800	97.550000
B	−170.570000	−179.200000	171.360000	−55.465000	5.695100	12.401000	106.606000
C	−201.660000	−61.000000	171.960000	−58.130000	6.230670	18.306000	141.354000
D	−239.860000	−247.760000	300.330000	−66.877000	6.885390	−12.233160	111.677000
E	−166.880000	−142.200000	153.160000	−55.066000	5.257900	−0.862400	140.425000
F	−210.970000	−88.990000	159.590000	−62.493000	7.755942	0.727800	119.188000
G	−225.880000	−114.140000	217.530000	−68.132000	5.880370	13.594400	103.596000
H	−216.380000	−33.040000	110.010000	−57.772000	4.046390	−0.383400	76.459000
I	−256.240000	−176.750000	262.440000	−71.867000	7.368280	−4.637400	102.291000

---

Molecular recognition

From the nine hit compounds, only seven compounds (A, B, D, F, G, H, I) are filtered through the PSA and these compounds are analyzed for molecular recognition with HTLV-1 PR crystal structure contours and new lead pharmacophore. These libraries of purvalanol-A are designed with respect to contours of receptor binding site. This molecular recognition fixes where the compound is correctly suitable for molecular environment. Understanding the pharmacophoric features and evaluate specific drug–receptor interactions with respect to contours provides strong binding of protein–ligand interactions. The interaction identifies important specific drug–receptor interactions between the purvalanol-A with benzyloxycarbonyl in the R⁴ position as having the pharmacophore features of acceptor – 6, donor – 3, hydrophobic – 3, aromatic – 4, which fits the best interaction with HTLV-1 PR; in other words, it fits best to the contours of HTLV-1 PR. Availability of benzyl group, hydrophobicity, and stereochemistry of certain functional groups are found to be important for inhibiting HTLV-1 PR, the molecular recognition of the ligand’s pharmacophore, which is the active part of the drug molecule, searches for contours in the receptor active site and binds to it with hydrogen bond interactions. Fig. 5 clearly shows HBAR, HBDR, HPR, and HPBR charged environments created by OPLS-2005 that occupy the benzyloxycarbonyl in R⁴ position having pharmacophoric feature of A – 6, D – 3, H – 3, R – 4, and so it is shown to have good docking score and binding energy. The R⁴ – pivaloyl does not fill the molecular recognition due to non-relation of binding contour and ligand pharmacophore sites (which is represented in Fig. 5); due to this non-relation, the R⁴ pivaloyl has lost its capacity to dock with HTLV-1 PR and has insufficient binding energy. So the molecular recognition is important in receptor–ligand interactions; in other words, recognition of binding contours of HTLV-1 PR recognize the perfect pharmacophores of new leads and gives better docking score and binding energy and compounds that passed through docking and binding energy analyzed with structure and ligand based descriptor analysis for theoretical analysis of activity with respect to contours have the ability to inhibit HTLV-1 PR, which is a target for anti-cancer drug design.

Fig. 5 Molecular recognization of protein contour interactions with respective ligand pharmacophoric sites ensures strong interactions and binding potential towards new leads of HTLV-1 PR inhibitors. (A) Surface contours present in HTLV-1 PR active site. (B) Pharmacophoric sites present in purvalanol-A. (C) Red color mesh protein contours interactions with hydrogen bond acceptors of ligands. (D) Blue color mesh protein contours interactions with hydrogen bond donors of ligands. (E) Hydrophilic region – green mesh interactions with hydrophilic and aromatic rings of ligands. (F) Hydrophobic region of yellow mesh in proteins interacts with hydrophobic regions of ligand pharmacophoric sites.

Molecular dynamics simulation

The ligand–receptor complex was further refined by MD simulations for 20 ns. The essential dynamic behavior of HTLV-1 PR protein in a water model was regulated up to 20 ns. The Root Mean Square Deviation (RMSD) of HTLV-1 PR backbone structure with respect to the initial conformation was calculated as a function of time period to assess the conformational stability of the protein during the simulations (Fig. 6). RMSD of the apo-form was ∼0.13 nm after 3.1 ns simulation and remained stabilized until the end of the simulation. The initial and final confirmation of the protein structure has been morphed with Chimera and analyzed. The backbone structure of the protein remains the same until the end of the simulation and the deviations are much less than 0.3 Å. The stable conformation obtained from dynamics studies will enhance the success rate of docking interactions and so we chose the average stable conformation for the computational part. The RMSD of apo (receptor alone) and the seven holo forms (receptor with ligands) during the MD simulation are shown in Fig. 6. Similarly, the average RMSD of the fluroenylmethoxycarbonyl (R¹ – A), phosphate (R² – B), benzyloxycarbonyl (R⁴ – D), sulfate (R³ – F), phenyloxycarbonyl (R⁴ – G), trifluroacetyl (R³ – H) and fluroenylmethoxycarbonyl (R⁴ – I) (seven holo-forms) showed ∼0.34 nm, ∼0.27 nm, ∼0.29 nm, ∼0.33 nm, ∼0.33 nm, ∼0.38 nm and ∼0.31 nm, respectively. The average RMSD of holo forms is higher than the apo form of HTLV-1 PR and both the bound and unbound forms of simulated structure remained same until the end of simulation. Ligand bound structures showed deviations, which shows that the new ligands are active inside the binding pocket. The hydrogen bonds formed are uniformly strong and hold the new lead compounds inside the HTLV-1 PR binding pocket (Fig. 7). The seven holo forms (receptor with ligands) during the MD simulation show strong binding and their average hydrogen bond formation is statistically calculated. The average H-bonds of fluroenylmethoxycarbonyl (R¹ – A), phosphate (R² – B), benzyloxycarbonyl (R⁴ – D), sulfate (R³ – F), phenyloxycarbonyl (R⁴ – G), trifluroacetyl (R³ – H) and fluroenylmethoxycarbonyl (R⁴ – I) are 3.8, 3.2, 3.3, 4.8, 2.8, 3.0 and 2.7 respectively. All the new lead compounds are show average H-bonds of ±3, which is higher than the known inhibitor purvalanol-A. These new compounds upon binding with stable potential energy support that new leads have strong potential to be held inside the binding pocket (Fig. 8). Potential energy fluctuations are not seen for any of the new compounds. On the whole, the new lead compounds are stable, show strong interactions, are very active and show binding potential towards HTLV-1 PR. These results suggest that binding of the ligand to the protein showed deviation from their initial position because of adjustments in their configuration but still remains to be bound within the active site of the protein.

Fig. 6 RMSD graph for HTLV PR apo and ligand complexes for the timescale event of 20 ns, shown with average mean variations.

Fig. 7 Hydrogen bond interactions of ligand-bound complex structures on the timescale of 20 ns and with average H-bond interactions.

Fig. 8 Potential energy of ligands exhibited during the molecular dynamics simulation for the timescale event of 20 ns.

Conclusion

Binding of small-molecule ligands to protein active sites is a key objective of drug design; theoretical interactions associated with receptors capture the effects of the ligand interaction with the protein active site with molecular recognition. The inhibition pattern of active sites are predicted with HTLV-1 PR and designed leads compounds of better interaction than purvalanol-A, which show up to −10.208523 of docking score, where purvalanol-A gives −6.702560. The molecular interaction reveals that SER55, GLY58, and ASN97 can also function as active sites from experimentally reported active sites as analyzes of both theoretical and experimental studies of active sites are important for drug design. Compounds of purvalanol-A attached R¹ – fluroenylmethoxycarbonyl, R² – phosphito, R³ – sulfato, R³ – trifluroacetyl, R⁴ – benzyloxycarbonyl, R⁴ – phenyloxycarbonyl, R⁴ – fluroenylmethoxycarbonyl show good molecular level interaction, which results in the ranking of top poses through binding energy, docking analysis, structure and ligand based descriptors and molecular dynamics simulations. The compounds that passed in docking interaction given as the best compounds may fail in binding energy, this is because of not considering the molecular level environment. Finally, the nine best leads of both docking interaction and binding energy were obtained based on the proper knowledge of both experimental and theoretical active sites, inhibiting HTLV-1 PR with designed small molecules based on requirement of contours present in active site. The inhibition of protein active sites is not just protein–ligand interaction, the actual phenomenon behind these studies involve with HBAR – A, HBDR – D, HPBR – R, HPR – H, and i.e. protein–ligand interaction deals with back end of contour of active site and ligand pharmacophore interactions. By these results, the interaction between the ligand and protein lies with contours and pharmacophore interactions. The ligand pharmacophore searches for an appropriate place in the binding site and it places well according to respective receptor contour surfaces. This molecular level recognition based on the ligand pharmacophore relies on respective contours, and finally this molecular attachment of ligand and receptors is fused with hydrogen bond formation.

Declaration of interest

None.

Acknowledgements

The research work is funded by CSIR and Sanjeev Kumar Singh thankfully acknowledges the CSIR for research funding (Ref. No. 37(1491)/11/EMR-II). Chandrabose Selvaraj gratefully acknowledges CSIR for a Senior Research Fellowship (SRF). Poonam Singh and Sanjeev Kumar Singh acknowledge CSIR-CDRI and Alagappa University for signing the Memorandum of Understanding. The CSIR-CDRI communication number allotted to this manuscript is 8816. Poonam Singh and Ankur Omer thankfully acknowledges CSIR-CDRI Director for his continuous support in research activities.

References

1. F. A. Proietti, A. B. Carneiro-Proietti, B. C. Catalan-Soares and E. L. Murphy, Oncogene, 2005, 24, 6058–6068.
2. P. Hollsberg and D. A. Hafler, N. Engl. J. Med., 1993, 328, 1173–1182.
3. R. F. Jarrett, J. Pathol., 2006, 208, 176–186.
4. K. Verdonck, E. Gonzalez, S. Van Dooren, A. M. Vandamme, G. Vanham and E. Gotuzzo, Lancet Infect. Dis., 2007, 7, 266–281.
5. Y. Satou and M. Matsuoka, J. Clin. Exp. Hematop., 2010, 50, 1–8.
6. C. Y. Hu, M. T. Lin, Y. C. Yang, J. L. Tang, L. H. Tseng, C. H. Wang, Y. C. Chen and C. S. Yang, J. Formosan Med. Assoc., 1998, 97, 101–105.
7. S. H. Nam, M. Kidokoro, H. Shida and M. Hatanaka, J. Virol., 1988, 62, 3718–3728.
8. P. Rucker, A. H. Horn, H. Meiselbach and H. Sticht, J. Mol. Model., 2011, 17, 2693–2705.
9. J. Tozser and I. T. Weber, Curr. Pharm. Des., 2007, 13, 1285–1294.
10. P. Boross, P. Bagossi, I. T. Weber and J. Tozser, Infect. Disord.: Drug Targets, 2009, 9, 159–171.
11. P. Singh, S. K. Singh, C. Selvaraj and R. K. Singh, Conversation 18, J. Biomol. Struct. Dyn., 2013, 31(1), 127.
12. M. Li, G. S. Laco, M. Jaskolski, J. Rozycki, J. Alexandratos, A. Wlodawer and A. Gustchina, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 18332–18337.
13. J. Kadas, P. Boross, I. T. Weber, P. Bagossi, K. Matuz and J. Tozser, Biochem. J., 2008, 416, 357–364.
14. J. Kadas, I. T. Weber, P. Bagossi, G. Miklossy, P. Boross, S. Oroszlan and J. Tozser, J. Biol. Chem., 2004, 279, 27148–27157.
15. C. Li, X. Li and W. Lu, Biopolymers, 2010, 94, 487–494.
16. S. B. Shuker, V. L. Mariani, B. E. Herger and K. J. Dennison, Chem. Biol., 2003, 10, 373–380.
17. L. Chen, J. K. Morrow, H. T. Tran, S. S. Phatak, L. Du-Cuny and S. Zhang, Curr. Pharm. Des., 2012, 18, 1217–1239.
18. S. Perot, O. Sperandio, M. A. Miteva, A. C. Camproux and B. O. Villoutreix, Drug Discovery Today, 2010, 15, 656–667.
19. N. K. Salam, R. Nuti and W. Sherman, J. Chem. Inf. Model., 2009, 49, 2356–2368.
20. H. J. Woo and B. Roux, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6825–6830.
21. D. L. Mobley, A. P. Graves, J. D. Chodera, A. C. McReynolds, B. K. Shoichet and K. A. Dill, J. Mol. Biol., 2007, 371, 1118–1134.
22. E. Agbottah, W. I. Yeh, R. Berro, Z. Klase, C. Pedati, K. Kehn-Hall, W. Wu and F. Kashanchi, AIDS Res. Ther., 2008, 5, 12.
23. SiteMap. Version 2.7: Schrödinger, LLC, New York, NY, 2013.
24. C. Liao, J. E. Park, J. K. Bang, M. C. Nicklaus and K. S. Lee, ACS Med. Chem. Lett., 2010, 1, 110–114.
25. T. Naumann and H. Matter, J. Med. Chem., 2002, 45, 2366–2378.
26. S. L. Dixon, A. M. Smondyrev and S. N. Rao, Chem. Biol. Drug Des., 2006, 67, 370–372.
27. D. Liszekova, M. Polakovicova, M. Beno and R. Farkas, PLoS One, 2009, 4, e6001.
28. C. Selvaraj, S. K. Singh, S. K. Tripathi, K. K. Reddy and M. Rama, Med. Chem. Res., 2012, 21, 4060–4068.
29. P. Vijayalakshmi, C. Selvaraj, S. K. Singh, J. Nisha, K. Saipriya and P. Daisy, J. Biomol. Struct. Dyn., 2013, 31, 561–571.
30. A. Volkamer, COMPASITES-Computer-aided active site analysis of protein structures.
31. D. C. Thompson, R. A. Denny, R. Nilakantan, C. Humblet, D. Joseph-McCarthy and E. Feyfant, J. Comput.-Aided Mol. Des., 2008, 22, 761–772.
32. Z. Li, H. Wan, Y. Shi and P. Ouyang, J. Chem. Inf. Comput. Sci., 2004, 44, 1886–1890.
33. Glide, version 6.1, Schrödinger, LLC, New York, NY, 2013.
34. W. Deng and C. L. Verlinde, J. Chem. Inf. Model., 2008, 48, 2010–2020.
35. R. A. Friesner, J. L. Banks, R. B. Murphy, T. A. Halgren, J. J. Klicic, D. T. Mainz, M. P. Repasky, E. H. Knoll, M. Shelley, J. K. Perry, D. E. Shaw, P. Francis and P. S. Shenkin, J. Med. Chem., 2004, 47, 1739–1749.
36. T. A. Halgren, R. B. Murphy, R. A. Friesner, H. S. Beard, L. L. Frye, W. T. Pollard and J. L. Banks, J. Med. Chem., 2004, 47, 1750–1759.
37. R. M. Shafreen, C. Selvaraj, S. K. Singh and S. K. Pandian, J. Mol. Recognit., 2013, 26, 276–285.
38. R. Farid, T. Day, R. A. Friesner and R. A. Pearlstein, Bioorg. Med. Chem., 2006, 14, 3160–3173.
39. C. Selvaraj and S. K. Singh, J. Biomol. Struct. Dyn., 2014, 32, 1333–1349.
40. C. K.-M. Chen, M. P. Hudock, Y. Zhang, R. T. Guo, R. Cao, J. H. No, P. H. Liang, T. P. Ko, T. H. Chang, S. C. Chang, Y. Song, J. Axelson, A. Kumar, A. H. Wang and E. Oldfield, J. Med. Chem., 2008, 51, 5594–5607.
41. H. S. Lee, J. W. Choi and S. J. Yoon, Genome inf., 2007, 5(1), 24–29.
42. P. K. Naik, M. Srivastava, P. Bajaj, S. Jain, A. Dubey, P. Ranjan, R. Kumar and H. Singh, J. Mol. Model., 2011, 17, 333–357.
43. A. Rizzi and A. Fioni, J. Chem. Inf. Model., 2008, 48, 1686–1692.
44. S. L. Dixon, A. M. Smondyrev, E. H. Knoll, S. N. Rao, D. E. Shaw and R. A. Friesner, J. Comput.-Aided Mol. Des., 2006, 20, 647–671.
45. N. Moitessier, C. Henry, B. Maigret and Y. Chapleur, J. Med. Chem., 2004, 47, 4178–4187.
46. J. R. Arnold, K. W. Burdick, S. C. Pegg, S. Toba, M. L. Lamb and I. D. Kuntz, J. Chem. Inf. Comput. Sci., 2004, 44, 2190–2198.
47. S. Pronk, S. Pall, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854.
48. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. Berendsen, GROMACS: fast, flexible, and free, J. Comput. Chem., 2005, 26(16), 1701–1718.
49. C. Selvaraj, J. Sivakamavalli, V. Baskaralingam and S. K. Singh, J. Recept. Signal Transduction Res., 2014, 34, 221–232.
50. C. Selvaraj, P. Singh and S. K. Singh, Appl. Biochem. Biotechnol., 2014, 172, 1790–1806.
51. A. W. Schuttelkopf and D. M. van Aalten, Acta Crystallogr., Sect. D: Biol. Crystallogr., 2004, 60, 1355–1363.
52. R. M. Beema Shafreen, C. Selvaraj, S. K. Singh and S. Karutha Pandian, J. Mol. Recognit., 2014, 27, 106–116.
53. C. Selvaraj, P. Singh and S. K. Singh, J. Mol. Recognit., 2014, 27, 696–706.
54. C. Selvaraj, P. Singh and S. K. Singh, J. Recept. Signal Transduction DOI:10.3109/10799893.2014.898659.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4mb00486h

---