Comparative study of the unbinding process of some HTLV-1 protease inhibitors using unbiased molecular dynamics simulations.

The HTLV-1 protease is one of the major antiviral targets to overwhelm this virus. Several research groups have developed protease inhibitors, but none has been successful. In this regard, developing new HTLV-1 protease inhibitors to fix the defects in previous inhibitors may overcome the lack of curative treatment for this oncovirus. Thus, we decided to study the unbinding pathways of the most potent (compound 10, PDB ID 4YDF, Ki = 15 nM) and one of the weakest (compound 9, PDB ID 4YDG, Ki = 7900 nM) protease inhibitors, which are very structurally similar. We conducted 12 successful short and long simulations (totaling 14.8 μs) to unbind the compounds from two monoprotonated (mp) forms of protease using the Supervised Molecular Dynamics (SuMD) without applying any biasing force. The results revealed that Asp32 or Asp32' in the two forms of mp state similarly exert powerful effects on maintaining both potent and weak inhibitors in the binding pocket of HTLV-1 protease. In the potent inhibitor's unbinding process, His66' was a great supporter that was absent in the weak inhibitor's unbinding pathway. In contrast, in the weak inhibitor's unbinding process, Trp98/Trp98' by pi-pi stacking interactions were unfavorable for the stability of the inhibitor in the binding site. In our opinion, these results will assist in designing more potent and effective inhibitors for the HTLV-1 protease.

Introduction

Human T-cell leukemia virus type 1 (HTLV-1) was discovered in 1980 as the first oncogenic retrovirus in the project "War on Cancer" in the United States [1]. According to the latest information, 5–10 million people are infected with this virus worldwide, and only 0.25–5% of them are affected by Adult T-cell Leukemia/Lymphoma (ATLL) and HTLV-1-associated myelopathy/tropical spastic paraparesis (HAM/TSP) [2], and also HTLV-1-associated ocular diseases. These diseases are known as HTLV-1 uveitis (HU) and ATL-related ocular [3]. Indeed, the reported numbers are not terrible, but there is no standard treatment for all types of diseases [4]. In addition, only a few regions were evaluated, and many unknown infected people could transmit the virus [5]. So, even low-risk areas are in danger because of Global Village. After HTLV-1 discovery, all its components were identified gradually, and its protease was discovered in 1989 [6]. HTLV-1 protease is a homodimer protein containing 125 residues in each subunit, which is one of the A2 family of aspartic proteases, with two critical aspartates in the catalytic dyad. This enzyme is essential for viral growth because it cleaves the Gag-Pro-Pol-Env polyprotein, a necessary viral replication component [7]. Since this part is vital for the viral life cycle, it is an interesting target for HTLV-1 demise.

Toward this end, many research groups in different countries succeeded in designing and synthesizing various compounds with inhibitory effects in the micromolar to nanomolar ranges [8, 9]. Finally, some German scientists considered the structural similarities between HTLV-1 and HTLV-3 (HIV) and determined the X-ray structure of Indinavir complexed with HTLV-1 protease, which is the only AIDS protease drug that has an inhibitory effect on HTLV-1 protease in low micromolar concentration. Unfortunately, this drug failed to be used to eradicate HTLV-1 [10]. After being frustrated with AIDS drugs, this team, in 2015, succeeded in synthesizing ten inhibitors that contain the most potent nonpeptidic inhibitor of HTLV-1 protease [11].

All reported HTLV-1 protease compounds only remain as inhibitors, and we do not have any specialized FDA-approved drug for inhibiting this virus. It is evident that experimental research alone is not sufficient. In silico methods, like unbiased molecular dynamics (UMD), are needed to provide valuable information for rational drug design, which is the primary goal of all researchers in this field. So, various research groups have used in silico methods for different studies, involving: carrying out the docking simulation to investigate the substrate-binding cavity of the IDH1 enzyme [12], useing dynamic simulations and binding free energy studies to design new potent inhibitors against CDK2 [13], prediction of the binding pocket in the interface of the aurora-A-TPX2 complex using molecular docking [14], finding specific naturally originated antagonist of the benzodiazepine binding site by performing docking experiments and molecular mechanics Poisson-Boltzmann surface area analysis [15]. MD simulation offers information about the reaction pathways of the ligand-protein complexes, and it has been considered by many research groups over these years and led to effective drug design [16, 17]. Therefore, besides the importance of one particular drug’s binding affinity to a target protein in the traditional drug design, the binding and unbinding processes and the residence time of the compound that interacts with the protein in each intermediate state are just as important. So by a complete understanding of the unbinding mechanism, we can uncover the key elements in the protein-ligand complex interactions [18], ligand flexibility, and solvation effects that are more critical in the rational drug design. The vital information will ultimately appear in a scenario with fully atomistic details [19]. For investigating unbinding pathways of inhibitors, some advanced MD simulation approaches like metadynamics and supervised metadynamics (suMetaD) simulation have been used before [20, 21], and one of the newest MD approaches is the supervised molecular dynamics (SuMD) method. This method performs simulation in replicas with fixed parameters. It gets information regarding metastable intermediate ligand-bound states using a tabu-like supervision algorithm. It is possible to fully unbind small molecules from their molecular targets without any biasing force or potential. In this regard, the SuMD has been utilized to discover the reaction pathways of various ligands in molecular targets [22].

The previous work of our team attempted to achieve Indinavir’s interaction with HIV protease and HTLV-1 protease as a comparison test to better understand the complexity of these key proteins [23]. In the following, we decided to examine the unbinding pathways of the most potent and one of the weakest HTLV-1 protease inhibitors retrieved from the last designed compounds using SuMD. In this study, the simulations will reveal the dynamic behavior of HTLV-1 protease as well as the strengths and weaknesses of the selected inhibitor. So, the results will guide the design of new potent HTLV-1 protease inhibitors. Also, this unbiased comparative study caused monitoring of all important factors in the two different inhibitor-protein interactions. So the outcomes were more reliable.

Methods

The X-ray crystallography structures of the HTLV-1 protease-ligand complex (PDB IDs: 4YDG, 4YDF [11]) were obtained from the Protein Data Bank. At first, for protein preparation, all protein missing residues and atoms in 4YDG were remodeled and fixed using UCSF Chimera software [24]. Then, for the ligand preparation, according to the practical information, the nitrogen of the pyrrolidine ring was protonated in both compounds [11] and parameterized by ACEPYPE using default settings (the GAFF atom type and BCC partial charges) [25]. HIV protease is an aspartic protease, and there were no studies on HTLV-1 protease, so we used a pattern from HIV research. In some research, protonation is not applied or not mentioned, or in some research, protonation state is applied in chain A. As we wanted a more complete outcome, we used protonation separately in both chains [26, 27]. So after preparing complexes, each catalytic Asp was considered separately as an ionization state based on the monoprotonated (mp) form of the catalytic dyad Asp32-Asp32′ in the active site [28]. Finally, we constructed our systems in GROMACS 2018 [29] using the OPLS all-atom force field [30] and with the TIP3P water model [31]. The considered holoproteins were centered in a triclinic box with a distance of 1 nm from each edge. The next step was to provide a 150 mM neutral physiological salt concentration, sodium, and chloride ions. Then all systems were relaxed in energy minimization using the steepest descent algorithm and reached Fmax of less than 1000 kJ.mol-1.nm-1. All covalent bonds by Linear Constraint Solver (LINCS) algorithm were constrained to maintain constant bond lengths [32]. The long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME) method [33]. The cut-off radii for Coulomb and Van der Waals (VdW) short-range interactions were set to 0.9 nm for all systems.

At last, the modified Berendsen (V-rescale) thermostat [34] and Parrinello-Rahman barostat [35] were applied for 100 and 300 ps, respectively, for keeping the system in stable environmental conditions (310 K, 1 Bar). To reach complete unbinds, we performed 12 separate replicas (three replicas for each type of mp form) with fixed duration times using the SuMD method with some modifications [36]. During the simulation, the distance between the center of masses of the ligands and selected residues was monitored until complete unbind occurred. This method is based on a tabu-like supervision algorithm without applying any human or non-human biasing force or potential. Herein, we set the center of mass (COM) of ligands as the first spot and the COM of the catalytic aspartic acids (Asp32, Asp32′) as the second spot and ran all simulations with a time window of 500 ps and a time step of 2 fs. After finishing each run, the frame with the longest distance between COMs was selected automatically to extend the next 500 ps simulation. These processes were continued until complete unbind was obtained, which is equal to a distance of 50 Å between the mentioned spots. Finally, all events in every concatenated trajectory file were investigated carefully with GROMACS utilities for data analysis. To picture the important interactions, we used UCSF chimera and used Daniel’s XL Toolbox (v7.3.4) to create plots [37], and using Matplotlib to show free energy landscapes (S1 Fig) [38]. The free energy landscapes plots were made base on three variables time, ligand RMSD, and protein RMSD. The ligand and protein RMSD values were selected because they were meaning full and had sharp changes as a function of time during unbindings. Analyzing these plots can reveal the stable states of inhibitors, as well as the residence time of inhibitors in each state over unbinding. Areas that tend to turn blue indicate that the inhibitor has been present in this area for a longer time.

Results and discussion

Since the only structural difference between compounds 9 and 10 is in the amino and nitro groups on the benzene ring (Fig 1A and 1B), compound 10 (Ki = 15 nM) is approximately 526 times more potent in complex with HTLV-1 protease [11]. Therefore, a proper understanding of the unbinding pathways of these compounds is vital to unveiling secrets that a minor structural difference can have a dramatic effect on inhibitory effects. Therefore, understanding the unbinding process of these inhibitors may provide an insight into the design of next-generation inhibitors.

Fig 1

The 2D structures of selected compounds were obtained from PDB.

A, Compound 9, the inhibitor in (PDB ID 4YDF). B, Compound 10, the inhibitor in (PDB ID 4YDG).

For a complete understanding of unbinding mechanisms of these two compounds, it is better to get more familiar with this less-known virus′s protease structure and inhibitors features at first. This homodimer protein has some particular regions with strategic effects in keeping ligands in the protein’s binding pocket that are obtained from the analysis of trajectories. The active site region (Leu31-Val39 and Leu31′-Val39′) contains catalytic dyad aspartate residues (Asp32, Asp32′) that are so important in protein-inhibitor interactions. The second essential region is the flaps (Val56-Thr63 and Val56′- Thr63′). The specific residues of Ala59-59′ consider as flap tips in the region of the flaps. Finally, Lateral Loops or 95S loops part of protease (Lys95-Gly102 and Lys95′- Gly102′) are other key regions in this aspartic protease (Fig 2A). For the inhibitors, both compounds have pi-pi intramolecular interactions. With more details, in compound 9, the nitrobenzene ring can form face-to-face pi-pi interaction with the benzene ring (Fig 2B), and in compound 10, the aniline ring can form T-shaped edge-to-face pi-pi interaction with the benzene ring (Fig 2C).

Fig 2

The 3D structure of HTLV-1 protease (PDB ID 4YDF) and intramolecular interactions of inhibitors.

A, All important domains of HTLV-1 protease: the green area is the active site region, the purple area is the lateral loops or 95S loops part of protease, the orange area is the flaps region, and the blue area is the flap tips part. B, Compound 9. C, Compound 10.

As mentioned, one of the essential parts of this protein is the region of the flaps, which showed high flexibility during simulations (Fig 3E–3H). So, during our simulations, four modes were observed for the flaps. Herein, we considered two factors to show these modes: the distance between COMs of Ala59 and Ala59′ (d1) and the second one is the distance between COMs of Ala59′ and Asp32′ (d2). The second factor can be even between COMs of Ala59 and Asp 32 due to flaps′ handedness opening. In the close form, the maximum amount of d1 and d2 are 10 and 15 Å, respectively (Fig 3A). In the semi-open form, the maximum amount of d1 and d2 is 14 and 20 Å, respectively (Fig 3B). In the open state, the minimum amount of d1 is 14 Å, and the maximum amounts of d1 and d2 are 20 Å (Fig 3C). In the wide-open form, d1 and d2 must be more than 20 Å (Fig 3D).

Fig 3

Different modes of the flap.

A, Close form of the flap. B, Semi-open form of the flap. C, Open form of the flap. D, Wide-open form of the flap. E, F, RMSF values of HTLV-1 protease in the 4YDG PDB code, during simulations. G, H, RMSF values of HTLV-1 protease in the 4YDF PDB code during simulations.

For our purpose, we had various simulations, that in total, we could have 12 successful unbindings (a total of 14.800 μs) for the two compounds from a minimum of 94 ns to a maximum of 4.4 μs in both mp forms. At last, for providing comprehensive information, all the events in each frame of trajectories were investigated carefully, and different analyses were performed on them.

Like the weak compound, we had two mp forms of the potent compound (AspH32 and AspH32′). In this regard, in the duration times of 4.4 μs (Fig 4A) and 260 and 305 ns (Fig 4B), which were in the chain A, Asp32 protonated state, we saw a uniform mechanism to unbind with some essential differences that caused a significant difference in one of simulation time. So, in the first state of rep1, 2, and 3 (Fig 4C–4E), Asp32′, which had salt bridge interaction with the positive charge of the pyrrolidine ring, plays a crucial role in preserving ligand in the binding pocket of protease. This acidic residue is essential because it is located almost at the bottom and center of the binding pocket. This residue considers a strategic residue due to the positive charge of pyrrolidine. Parallel to that, His66′ by cation–pi interaction [39, 40] with an aniline ring was the second important preserving factor. In addition, Ala59′ in the flap tip by forming H-bond with the atom of O10 (Fig 4H) and also Ala35′ by VdW interactions with a benzene ring, Asp36′ by forming H-bond with an aniline fragment in the active site, and finally Ile100′ in 95S loop (Fig 4J–4L), with VdW interaction, blocked all the exit routes, like the fence (Fig 4F). As mentioned before, His66′ was the second essential residue in this state, which was a supporter of Asp32′ to fix the inhibitor in the binding pocket. As time passing, Asp32′ loosed its superpower of preserving, and the inhibitor entered the second intermediate state. In this state, Lys95′, forming a hydrogen bond with the atom of O10, along with His66′ cation-pi interaction with the benzene ring, was a third essential residue. This residue increased protein-ligand interactions time in the rep1 and was absent in rep2 and rep3 (Fig 4G). According to significant differences in replicas simulation times, the effects of the Lys95′ hydrogen bond (Fig 4I) appear more pronounced. Finally, ligand pi-pi intramolecular interactions, observed during the whole time of simulations (Fig 5A–5C), slowly weakened all-important protein-inhibitor interactions. The critical point was that, over the entire simulation time, flaps positioning impacted the ligand’s behaviors, so the exit process started when the flaps began to open, and Ala59′ loosed its effect (Fig 5D–5F) in the second intermediate state gradually with the help of water mediation (Fig 6A–6C).

Fig 4

The details of compound 10 unbinding pathways in complex with HTLV-1 protease when Asp32 of chain A was protonated in three replicas.

A, RMSD value of the ligand from binding pose to complete unbinding in the rep1. B, RMSD values of the ligand from binding pose to complete unbinding in the rep2 and rep3. C, D, and E, The free energy landscape of rep1, 2, and 3 during the unbinding process (state (S), intermediate state (I), unbound (U)), respectively, which was calculated by using "gmx sham". F, The interactions between the ligand and essential residues in the binding pose of rep1, 2, and 3. G, The new interactions between the inhibitor and particular residues in the second intermediate state of rep1. H, Hydrogen bond numbers of Asp36′ and Ala59′ with the inhibitor in rep1, 2, and 3. I, Hydrogen bond numbers of Lys95′ with the inhibitor in rep1. J, K, and L, The average of most important interaction energies of the protein-ligand complex in rep1, 2, and 3, respectively.

Fig 5

The details of distances between particular parts of compound 10 complexed with HTLV-1 protease when Asp32 of chain A was protonated in three replicas.

A, B, and C, The distance between COMs of both aniline rings and benzene rings, which were in a position that could form pi-pi intramolecular interactions in all replicas. These plots prove that these fragments were so close together during the simulation. D, E, and F, The distance between COMs of Ala59 and Ala59′, and also Asp32′ and Ala59′ in all replicas (these plots should be checked along with Fig 3).

Fig 6

The interaction energies plots of compound 10 in complex with HTLV-1 protease when Asp32 of chain A was protonated in three replicas.

A, B, and C, The total VdW and electrostatic interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

Conversely, in the other state of protonation (AspH32′), we saw a uniform pathway that was dissimilar to the previous unbinding proccess with the different lengths of times involving: 94, 320, and 790 ns (Fig 7A and 7B). In the first state (S1 Fig) of these pathways (Fig 7C–7E), Asp32 was as important as expected. Asp36 and Asp36′, Leu57, and Ala59′ are the residues that acted as auxiliary agents (Fig 7J–7L) to the pivotal amino acid (Asp32). At the first simulation times, along with the salt bridge of Asp32 and pyrrolidine fragment (Fig 7F), both aniline rings had H-bonds with Asp36 and Asp36′ in the active site (Fig 7H). Along with these residues, Leu57 and Ala59′ formed a hydrogen bond with an aniline fragment and O41 atom of inhibitor, respectively (Fig 7G and 7I). In the following, in the lack of His66 and Lys95 effects, after time passing with the help of pi-pi ligand intramolecular interactions and water molecules effect (Fig 8A–8C), active site and flaps’ important residues lost their effects, and full unbind was observed between the flaps (Figs 8D–8F and 9A–9C).

Fig 7

The details of compound 10 unbinding pathways in complex with HTLV-1 protease when Asp32′ of chain B was protonated in three replicas.

A, RMSD values of the ligand from binding pose to complete unbinding in the rep1 and rep2. B, RMSD value of the ligand from binding pose to complete unbinding in the rep3. C, D, and E, The free energy landscape of rep1, 2, and 3 during the unbinding process (state (S), intermediate state (I), unbound (U)), respectively, which was calculated by using "gmx sham". F, The interactions between the ligand and essential residues of the active site in the binding pose of rep1, 2, and 3 G, The interactions between the ligand and essential residues of the region of the flaps in the binding pose of rep1, 2, and 3. H, Hydrogen bond numbers of Asp36 and Asp36′ with aniline fragment in rep1, 2, and 3. I, Hydrogen bond numbers of Leu57 and Ala59′ with the inhibitor in rep1, 2, and 3. J, K and L, The average of most important interaction energies of the protein-ligand complex in rep1, 2, and 3, respectively.

Fig 8

The details of distances between particular parts in compound 10 in complex with HTLV-1 protease, when Asp32 of chain B was protonated, in three replicas.

A, B, and C, The distance between COMs of both aniline rings and benzene rings, which were in a position that could form pi-pi intramolecular interactions in all replicas. These plots prove that these fragments were so close together during the simulation. D, E, and F, The distance between COMs of Ala59 andAla59′ and also Asp32′ Ala59′ in all replicas (these plots should be checked with Fig 3).

Fig 9

The interaction energies plots of compound 10 in complex with HTLV-1 protease when Asp32 of chain B was protonated in three replicas.

A, B, and C, The total VdW and electrostatic interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

Depending on close and open flaps and mp states, we had different unbinding mechanisms for the weaker inhibitor. Accordingly, compound 9 was unbound in 148 ns (Fig 10A), 3.5 μs, and 3 μs (Fig 10B) when Asp 32 of chain A was protonated. In the first state (S1 Fig) (Fig 10C) of the rapid unbinding pathway (rep1), a repulsive force occurred between the pyrrolidine ring of ligand and AspH32 of the binding pocket, and because of no attractive interactions in this area, AspH32 forced the ligand to push out. In addition to this interaction, VdW interaction between both nitrobenzene rings and one of the benzene rings of inhibitor and Leu57, Gly58, and Ala59 in the close flap region and pi-pi stacking interaction of Trp98′ and nitrobenzene fragment and also pi-alkyl interaction of Ile100′ with the benzene ring, were other protein-inhibitor significant interactions, which were not potent enough to prevent from repulsive interaction effect (Fig 10F and 10K). In the two other longer simulations, compound 9 was unbound in 3.5 μs in wide-open flaps (rep2) and 3 μs in close and semi-open flaps (rep3).

Fig 10

The details of compound 9 unbinding pathways in complex with HTLV-1 protease when Asp32 of chain A was protonated in three replicas.

A, RMSD value of the ligand from binding pose to complete unbinding in the rep1. B, RMSD values of the ligand from binding pose to complete unbinding in the rep2 and rep3. C, D, and E, The free energy landscape of rep1, 2, and 3 during the unbinding process (state (S), intermediate state (I), unbound (U)), respectively, which was calculated by using "gmx sham". F, The interactions between the ligand and essential active site residues in the rep1. G, The interactions between the ligand and essential active site residues in the rep2. H, The interactions between the ligand and essential active site residues in the binding pose of rep3. I, The new interactions between the inhibitor and particular residues in the second intermediate state of rep3. J, The distance between COMs of pyrrolidine ring and Asp36 and Asp32′ in rep3, to show after 2us of simulation this fragment get closer to Asp36 and get farther from Asp32′. K, L, and M, The average of most important interaction energies of the protein-ligand complex in rep1, 2, and 3, respectively.

Similarly, Asp32′ was the most important amino acid with its salt bridge and the only common point in both pathways. In the first state of rep2 (S1 Fig) (Fig 10D), due to handedness opening, only one of the flaps had forward and backward motions, so Leu57′, Gly58′, Ala59′ by VdW interactions kept the ligand in exposing to Asp32′. Also, in this state, Trp98 in the lateral loop built up pi-pi stacking interaction with the nitrobenzene ring of the ligand, and Trp98′ built up pi-pi stacking interaction with the benzene ring of another side of the inhibitor (Fig 10G and 10L). So even with enough space for the exit, the inhibitor was still in blockage. These important protein-inhibitor interactions were maintained until the effect of the Asp32′ became faded, and other agents lost their effect. Unexpectedly, the interesting point was that the complete unbinding process does not occur from the region of the flaps. In the rep3 pathway, that the flaps were close or semi-open the whole time, from the first state (Fig 10E), not only Asp32′ was necessary, and Asp36 in a close position to Asp32′ was powerful too (Fig 10M).

On the other hand, during the first two states, Asp36, by forming pi-anion interaction [41] with nitrobenzene fragment, was momentous as a second ligand preserving residue (Fig 10H), which was promoted to the first important factor in the following intermediate state by replacing pi-anion interaction with the salt bridge with pyrrolidine fragment (Fig 10I). From a holistic view, even though Asp32′ was more critical for protein, it was effective until the second intermediate state or until 2 μs, but Asp36 (Fig 10J) was effective until complete unbind. Actually, the ligand in all replicas showed face-to-face pi-pi intramolecular interactions between mentioned fragments that caused weakened important protein-ligand interactions gradually with the help of the water mediation effect (Fig 11A–11C). Finally, for the flaps behaviors in all replicas, we saw a new opening form for the rep2 as it was opened from chain A (Fig 11E), and for rep1 and rep3 wide opening (Fig 11D and 11F) was not seen until complete unbound (Fig 12A–12C).

Fig 11

The details of distances between particular parts in compound 9 in complex with HTLV-1 protease, when Asp32 of chain B was protonated, in three replicas.

A, B, and C, The distance between COMs of both nitrobenzene rings and benzene rings, which were in a position that could form pi-pi intramolecular interactions in all replicas. These plots prove that these fragments were so close together during the simulation. D, The distance between COMs of Ala59 andAla59′ and also Asp32′ Ala59′ in the rep1 (these plots should be checked with Fig 3). E, The distance between COMs of Ala59 andAla59′ and also AspH32 Ala59 in the rep2. F, The distance between COMs of Ala59 andAla59′ and also Asp32′ Ala59′ in the rep3. G and H, The total interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

Fig 12

The interaction energies plots of compound 9 in complex with HTLV-1 protease when Asp32 of chain A was protonated in three replicas.

A, B, and C, The total VdW and electrostatic interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

On the contrary, when the Asp32 of chain B was protonated, we saw the same mechanism during the 1.2 μs, 410, and 450 ns of simulations (Fig 13A and 13B). In the first state of these replicas (Fig 13C–13E), the ligand was surrounded by interactions of some residues in both chains (Fig 13G–13I). Asp32 had salt bridge interaction with pyrrolidine fragment as a most important interaction. In more detail, this fragment also had VdW interaction with Gly34, one of the nitrobenzene rings was in VdW interactions with Leu57′ and Gly58′ in the flaps regions for the other fragments. The benzene rings were in important interactions involving: pi-pi interaction with Trp98 and pi-alkyl interaction with Ile100 on one side, and pi-pi interactions with Trp98′ and pi-alkyl interaction with Ile100′ on the other side (Fig 13F). It may be due to the high number of important factors; it seems that compound 9 is potent, but except for Asp32 other agents did not have any significant effect. So, they could not keep the ligand after disappearing the Asp32 effect. Thus, as time passed, intramolecular interactions of ligand and water mediation contributed to full unbinding in these three replicas (Fig 14A–14C). Ultimately for the flaps effects, in rep1, the flaps showed high motions, and even though the flaps were wide open (Fig 14D), the full unbind did not occur from this region. In rep2 and 3, the inhibitor unbounded between semi-open flaps forms (Figs 14E, 14F and 15A–15C).

Fig 13

The details of compound 9 unbinding pathways in complex with HTLV-1 protease when Asp32 of chain B was protonated in three replicas.

A, RMSD value of the ligand from binding pose to complete unbinding in the rep1. B, RMSD values of the ligand from binding pose to complete unbinding in the rep2 and rep3. C, D, and E, The free energy landscape of rep1, 2, and 3 during the unbinding process (state (S), intermediate state (I), unbound (U)), respectively, was calculated using "gmx sham". F, The interactions between the ligand and essential active site residues in the rep1, 2, and 3. G, H, and I, The average of most important interaction energies of the protein-ligand complex in rep1, 2, and 3, respectively.

Fig 14

The details of distances between particular parts of compound 9 complexed with HTLV-1 protease when Asp32 of chain B was protonated in three replicas.

A, B, and C, The distance between COMs of both nitrobenzene rings and benzene rings, which were in a position that could form pi-pi intramolecular interactions in all replicas. These plots prove that these fragments were so close together during the simulation. D, E, and F, The distance between COMs of Ala59 andAla59′ and also Asp32′ Ala59′ in all replicas (these plots should be checked with Fig 3).

Fig 15

The interaction energies plots of compound 9 in complex with HTLV-1 protease when Asp32 of chain B was protonated in three replicas.

A, B, and C, The total VdW and electrostatic interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

Conclusion

The atomistic details of unbinding pathways of selected inhibitors in all replicas with various times, the importance of Asp32′ in chain A protonation state and Asp32 in chain B protonation state are pretty straightforward. Due to its strategic position, this effective residue could play a critical role in keeping the ligand in the binding pocket for a long time, so the more exposed to Asp32 or Asp32′, the more inhibitory effects. The pyrrolidine fragment was held well by Asp32 or Asp32′ from the native binding pose of the two compounds, which correlates with experimental research [11]. Thus the interactions of other fragments with other residues in different protein regions caused significant differences.

Herein, we cannot conclude which state of protonation occurs, so with our obtained information for the potent compound in chain A protonation state, His66′ with its cation–pi interaction with an aniline ring of inhibitor was a perfect supporter to Asp32′. This residue’s effect was absent in the other form of protonation state and caused a significant difference in simulation time. In the weak inhibitor unbinding pathways, Trp98 and Trp98′ with pi-pi interactions, due to their close position to one of the exit areas, were not good supporters for Asp32 or Asp32′, like His66′. His66′, due to its far position from the bottom and center of the binding pocket, could fix aniline fragments and decrease ligand fluctuations. The two mentioned tryptophan were closer to the essential aspartic acids, and there was enough space for ligand fluctuations. For this reason, Asp36 in the active site that was close to the exit area could be a competitor with Asp32′ and was not a good interaction for keeping the ligand in the binding pocket. Similarly, attenuating effect of Trp98/Trp98′ residues in unbinding pathways of the weak inhibitor correlates with our other research result [23]. These residues’ interactions are unfavorable to Indinivar’s stability in complex with HTLV-1 protease and result in it being a weak inhibitor. As we said before, both compounds had intramolecular interactions that caused weakening critical protein-ligand interactions as time passing. These two compounds did not have the same intramolecular interactions type, so in the weak inhibitor, face-to-face pi-pi interactions resulted in the loss of significant pi interactions with the protein. However, the potent inhibitor could have formed more important pi interactions with the protein along with intramolecular interactions. Overall, this obtained information is valuable for designing a new generation of inhibitors against this molecular target.

Because no similar simulation has been performed on the interaction of these inhibitors with HTLV-1 protease, we were only able to compare the results of native states with the crystallographic binding poses. So, except Asp32, Asp32′, in the potent inhibitor and Trp98, and in the weaker inhibitor, Asp36, Trp98 and Trp98′ were the only common critical residues [11].

The free energy landscape plots of all replicas.

(TIF)

Click here for additional data file.

(TIF)

Click here for additional data file.

15 Mar 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Sandipan Chakraborty

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Partly

Reviewer #3: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: N/A

Reviewer #3: No

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Comments;

1. Abstract is very poorly written. Rewrite the abstract according to findings.

2. Clearly define the aim and objectives of the study in the last paragraph of the introduction section.

4. In the Introduction section the author should refer to the research paper and comment on recent in-silico techniques. It will be good information for the readers. I would like to recommend several papers, among many others, providing further explanation on this topic:PMID: 27194485 PMID: 33749525 PMID: 32448055 PMID: 31980008 PMID: 33065246 PMID: 32447145 PMID:  34717229

4. Authors have advised redrawing all the interaction energy graphs in reverse order.

5. The interaction energy values appear to be much larger than those of typical small molecule inhibitors. I think the interaction energy is insufficient to determine the unbinding mode of the ligand accurately. Is there a possibility of overestimation such as incorrectness of reweighing or insufficient sampling of the unbound state?

6. Authors have provided insufficient data. Authors have advised to perform binding studies first before exploring the unbinding.

7. Provide a figure showing the unbinding along with the simulation time.

8. Provide pull parameters in the methods section. Also, briefly explain the general theory of Supervised Molecular Dynamics. How the molecules are pulled out of the pocket, in what direction were the ligands pulled, etc.

9. Define abbreviations in the Abstract section.

10. Why OPLS all-atom force field was used?

11. Authors have not provided the data of external pulling force and contacts to show the unbinding pathway of the selected molecules.

12. Results and discussion lacks data and needs to be elaborated in comparison with other computational data based on similar studies.

13. Authors have advised checking the quality of the minimized structure of the 1HTLV-1 and the authors should carry out additional docking studies with experimentally known inhibitors and compare the computational inhibition values with the experimental values. Without any experimental support or validation studies, the in silico binding free energy calculations may lead us totally wrong results, and the whole work may be nonsense.

14. Overall, the study is incomplete and requires more robust analyses to validate the findings experimentally or computationally. In addition, the manuscript is scientifically unsound and not suitable to publish in this journal without properly validated studies.

Few minor comments;

“We had two mp forms of the potent compound (AspH32 and AspH32′) like the weak compound. In this regard, in the duration times of 4.4us” Make corrections to the unit of time.

What is mp state? The authors have advised to provide abbreviations at their first use.

“In the weak inhibitor unbinding pathways, Trp98 and Trp98′ with pi-pi interactions, due to their close position to one of the exit areas were not good supporters for Asp32 or Asp32′, like His66′, His66′,” Author should correct the notation of His residue.

Authors have calculated the total interaction energy of compound 9 and 10 with htlv but did not mention anything about it in the text. Merely displaying in figures does not make sense.

“We had two mp forms of the potent compound (AspH32 and AspH32′) like the weak compound.” Should mention properly which is potent and each compound.

Authors already mentioned that one compound is better than other and they are just investigating the mechanism of unbinding. In MS authors tells about two mp form of compounds (potent and weak). Is the mp forms be taken into account for getting Ki values for these compounds?

In my view, the results obtained in this study are worthy for publication. The manuscript needs major essential revision before publication. I would like to overview the revised version of the manuscript before it accept for publication.

Reviewer #2: 1. What purpose do authors have selected for two different structures of HTLV-1 protease for this study?

2. Authors have used co-crystalized ligands in this study. What is the novelty of the proposed study?

3. From the literature it was reported that the amino acid Met37 played a key role in the binding of inhibitors in the active site of HTLV-1 protease. Did authors have noticed interaction with this amino acid in this study?

4. Authors have mentioned the uniform mechanism to unbind with some important differences: what are key features or differences noticed during the simulation period?

5. Authors have mentioned that compound shows 526 times more potent in complexes with HTLV-1 than compound 9, is the substitution of amino group of this compound enhances the activity of compound 10?

6. Does the amino group of the compound 10 show any interaction with the important amino acid residues of the HTLV-1 protease?

7. Authors have mentioned that the positive charge of the pyrrolidine ring plays a crucial role in preserving ligands. Both compound 9 and 10 have pyrrolidine rings in their structure, then how compound 9 shows weaker activity than compound 10? Does any specific mechanism play a role in activity of compound 10?

8. In the main text it was mentioned that the figure 4 c-e was the first state of rep1, 2, and 3, but in figure it was mentioned as ligand RMSD. Make it correct.

9. In figure 5b and c, both lines were depicted as Rep2 and Rep3-aniline-benzene. Then how does it show the difference in the distance? What are the main differences noticed in these structures?

10. Is the author using compound 9 as the weaker inhibitor in this study?

Reviewer #3: The manuscript is an in-depth study but needs extensive revisions for publications. Statistical analyses or any convergence tests are missing. It is well written, but figures need to make much better to look like as publishable format. Details review have been attached in the 'reviewer's feedback.docx' file.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: Reviewers feedback.docx

Click here for additional data file.

28 Apr 2022

Dear Prof. Sandipan Chakraborty;

I am pleased to mail you the revised article titled "Comparative study of the unbinding process of some HTLV-1 protease inhibitors using Unbiased Molecular Dynamics simulations". All comments have been answered upon reviewers' recommendations, and the corrections were made in manuscript # PONE-D-22-00753. It is our pleasure to hear your feedback and any other suggestions relevant to the paper.

Sincerely,

Hassan Aryapour

Reviewer #1:

Comment #1:

The abstract is very poorly written. Rewrite the abstract according to findings.

Answer #1:

The abstract section was improved and rewritten.

Comment #2:

Clearly define the aim and objectives of the study in the last paragraph of the introduction section.

Answer #2:

The requested information was added in the last paragraph of the introduction section (lines 100-104)

Comment #3:

In the Introduction section the author should refer to the research paper and comment on recent in-silico techniques. It will be good information for the readers. I would like to recommend several papers, among many others, providing further explanation on this topic: PMID: 27194485 PMID: 33749525 PMID: 32448055 PMID: 31980008 PMID: 33065246 PMID: 32447145 PMID: 34717229

Answer #3:

Some of your mentioned references were added as an example of some in silico methods (lines 72-79)

Comment #4:

Authors have advised redrawing all the interaction energy graphs in reverse order.

Answer #4:

In this study, we examined unbinding pathways of mentioned inhibitors, so at first, the inhibitors were in a binding position with the highest amount of interaction energies. During the simulation, the weakening of protein-ligand interactions occurred. The inhibitors started to leave the binding site and eventually full unbind occured when the total amount of protein-ligand interaction energies reached zero. So all the interactions energies plots' trend is correlated with the unbinding process from the bound state to the unbound state. This kind of analysis also used in the other studies:

https://doi.org/10.1371/journal.pone.0263251,https://doi.org/10.1371/journal.pone.0257916 https://doi.org/10.1371/journal.pone.0251910

Comment #5:

The interaction energy values appear to be much larger than those of typical small molecule inhibitors. I think the interaction energy is insufficient to determine the unbinding mode of the ligand accurately. Is there a possibility of overestimation such as incorrectness of reweighing or insufficient sampling of the unbound state?

Answer #5:

Total interaction energies plots are obtained from the sum of the coulomb and Lennard-jones contribution values, so the total energies are high. To better understand the complementary analysis, we calculate the average of the most important interaction energies (the amount of these plots is not much as the Total interaction energies plots).

Investigating the energy plots was not the only analysis for determining unbinding processes. We also used the free energy landscape plots obtained using "gmx sham", based on the ligand and protein RMSD values. Eventually, all preliminary information was investigated with trajectory files to advanced details obtain.

Comment #6:

The authors have provided insufficient data. Authors have advised to perform binding studies first before exploring the unbinding.

Answer #6:

Along with investigating the unbinding pathways of these two inhibitors we also investigated the binding pathways of these inhibitors that we are still writing the article. The information gathered in the binding process has covered the results of the unbinding process. Also there some valuable research regarding only unbinding pathways investigation such as: DOI: 10.1038/srep11539 , DOI: 10.1126/sciadv.1700014, DOI: 10.1073/pnas.1424461112.

Comment #7:

Provide a figure showing the unbinding along with the simulation time.

Answer #7:

In the Table of Contents graphic section, we provided a schematic picture of one of the compound 9 and 10 unbinding pathways as an example. RMSD values of the ligand from binding to complete unbinding positions are also available ( Fig4 a,b). RMSD plots show how ligand location changes with time as compared to their crystallographic binding poses.

Comment #8:

Provide pull parameters in the methods section. Also, briefly explain the general theory of Supervised Molecular Dynamics. How the molecules are pulled out of the pocket, in what direction were the ligands pulled, etc.

Answer #8:

We use a simulation method (SuMD) entirely different from “umbrella sampling” to pull ligands from their binding sites. A brief explanation of the SuMD method is available on lines 91-97. The directions of inhibitor unbinding are available in the Result and Discussion section, for example: " In rep2 and 3, the inhibitor unbounded between semi-open flaps forms (Figure 14e, 14f) (Figure 15a, 15b, 15c)".

Comment #9:

Define abbreviations in the Abstract section.

Answer #9:

The requested correction was done in the Abstract section (line 25).

Comment #10:

Why OPLS all-atom force field was used?

Answer #10:

Because we used the OPLS force field for binding simulations, the unbinding simulations were also done by it to reproduce the data. Also, in some other similar research involving: https://doi.org/10.1093/bioinformatics/btaa565, DOI: 10.1126/sciadv.1700014, this force field is used. In OPLS, surfaces (parameters) are fit to experimental data sets MM PES, so this FF is perfect for describing the properties of small molecules in their condensed states. OPLS-AA is entirely open-source, so it can be applied in many MD programs; particularly, the most popular MD program has already integrated the OPLS-AA into it.

Comment #11:

Authors have not provided the data of external pulling force and contacts to show the unbinding pathway of the selected molecules.

Answer #11:

The process of unbinding from binding pose to the full unbind is dependent on a set of factors, including aspects of protein-ligand interactions, water mediation, protein, and ligand fluctuation. In the Result and Discussion section, we explained the atomistic details of all of these factors with different analyses. As we answered in “Comment#8” we did not use pull parameters.

Comment #12:

Results and discussion lack data and need to be elaborated in comparison with other computational data based on similar studies.

Answer #12:

In the other computational research, the study's goal differed. Some of our data were new, so they could not be compared, except our other team's research that works on unbinding pathways of Indinavir in complex with HTLV-1 protease. Previously, we discussed our correlation of results in the conclusion section "These residues' interactions are unfavorable to Indinivar's stability in complex with HTLV-1 protease and result in it being a weak inhibitor".

Comment #13:

Authors have advised checking the quality of the minimized structure of the 1HTLV-1 and the authors should carry out additional docking studies with experimentally known inhibitors and compare the computational inhibition values with the experimental values. Without any experimental support or validation studies, the in silico binding free energy calculations may lead us totally wrong results, and the whole work may be nonsense.

Answer #13:

In the other project, we investigated the binding pathways of the two inhibitors, which will be available soon. By combining the results of this article and binding study results, we can have a more positive attitude towards drug design.

Only crystallographic binding pose information is analyzed in the experimental research, but we study the interactions during the pathways. As you run more simulations, you will discover more pathways. However, we added the results compared with experimental research in the conclusion section (lines 433-434 and 460-467) to better understand.

Few minor comments;

We had two mp forms of the potent compound (AspH32 and AspH32′) like the weak compound. In this regard, in the duration times of 4.4us" Make corrections to the unit of time.

What is mp state? The authors have advised to provide abbreviations at their first use.

"In the weak inhibitor unbinding pathways, Trp98 and Trp98′ with pi-pi interactions, due to their close position to one of the exit areas were not good supporters for Asp32 or Asp32′, like His66′, His66′," Author should correct the notation of His residue.

Authors have calculated the total interaction energy of compound 9 and 10 with htlv but did not mention anything about it in the text. Merely displaying in figures does not make sense.

"We had two mp forms of the potent compound (AspH32 and AspH32′) like the weak compound." Should mention properly which is potent and each compound.

Authors already mentioned that one compound is better than other and they are just investigating the mechanism of unbinding. In MS authors tells about two mp form of compounds (potent and weak). Is the mp forms be taken into account for getting Ki values for these compounds?

Answer:

A unit correction was done. The abbreviation of mp was added in the Abstract section. The correction was done in "In the weak inhibitor unbinding pathways, Trp98 and Trp98′ with pi-pi interactions, due to their close position to one of the exit areas were not good supporters for Asp32 or Asp32′, like His66′".

Total interaction energies plots are obtained from the sum of the coulomb and Lennard-jones contribution values. We have maximum interaction energies in the bound state of the protein-ligand complex. As time passes and protein-ligand interactions weaken, the values of the plot increase and eventually reach zero until full unbinding. So all the interactions energies plots' trend is correlated with the unbinding process from the bound state to the unbound state. So we used these plots to prove full unbind occurrence by energy values.

Aspartic proteases tend to cleave substrate by monoprotonation of aspartates in the catalytic dyad. Since HIV protease is an aspartic protease, and unfortunately, there was no typical study about HTLV-1 protease, we decided to take a pattern from HIV research. So in some research, they do not apply protonation or do not mention it, or in some research, they apply protonation state in chain A. Due to we wanted to have more complete outcome, we decided to use protonation in both chains separately.

We mentioned in the conclusion section, "Herein, we cannot conclude certainly which state of protonation actually occurs". In the presence of an inhibitor, the kind of protonation that will occur can only be determined when exposed to the active site. Because of this, we decided to discuss protonation more cautiously and, instead of judging, report the results.

Reviewer #2:

Comment #1:

What purpose do authors have selected for two different structures of HTLV-1 protease for this study?

Answer #1:

As we mentioned in the Result and Discussion section, "Since the only structural difference between compounds 9 and 10 is in the amino and nitro groups on the benzene ring (Figure 1a, 1b), compound 10 (Ki = 15 nM) is approximately 526 times more potent in complex with HTLV-1 protease. Therefore, a proper understanding of the unbinding pathways of these compounds is vital to unveiling secrets that a minor structural difference can have a dramatic effect on inhibitory effects." We could monitor critical factors from a comparative perspective, which made the results more reliable for designing new inhibitors. The related section improved in the manuscript.

Comment #2:

The authors have used co-crystallized ligands in this study. What is the novelty of the proposed study?

Answer #2:

There is no FDA-approved drug for this target, so we decided to work on designed inhibitors. In our view, it is critical to determine the negatives and the positives of the current potent inhibitor. We also compared the result with a weak inhibitor to achieve a complete insight for designing next-generation inhibitors.

Comment #3:

From the literature, it was reported that the amino acid Met37 played a key role in the binding of inhibitors in the active site of HTLV-1 protease. Did authors have noticed interaction with this amino acid in this study?

Answer #3:

In our replicas, Met37 effects were not significant as critical residues. Important interaction energies values are available in energy plots.

Comment #4:

Authors have mentioned the uniform mechanism to unbind with some important differences: what are key features or differences noticed during the simulation period?

Answer #4:

The Phrase "we saw a uniform mechanism to unbind with some important differences that caused a significant difference in one of simulation time" was used for compound 10 unbinding pathways in complex with HTLV-1 protease when Asp32 of chain A was protonated. We mentioned "In this state, Lys95′ by forming a hydrogen bond with the atom of O10, along with His66′ cation-pi interaction with the benzene ring, was a third essential residue" (lines 231-232) and "According to significant differences in replicas simulation times, the effects of the Lys95′ hydrogen bond (Figure 4i) appear more pronounced" (lines 235-236) as a key factor that caused significant different in simulation time in the Result and Discussion section.

Comment #5:

Authors have mentioned that the compound shows 526 times more potent in complexes with HTLV-1 than compound 9, is the substitution of the amino group of this compound enhances the activity of compound 10?

Answer #5:

Yes, it is. The inhibitor 10 could be able to interact with His66′ and Lys95′, which are essential.

Comment #6:

Does the amino group of the compound 10 show any interaction with the important amino acid residues of the HTLV-1 protease?

Answer #6:

Yes, it does. Asp36′ by forming H-bond when Asp32 of chain A was protonated, and Asp36 and Asp36′ by forming H-bonds when Asp32 of chain B was protonated with the amino group, which is available in the text line 226 and line 270-274.

Comment #7:

Authors have mentioned that the positive charge of the pyrrolidine ring plays a crucial role in preserving ligands. Both compound 9 and 10 have pyrrolidine rings in their structure, then how compound 9 shows weaker activity than compound 10? Does any specific mechanism play a role in activity of compound 10?

Answer #7:

We mentioned in the Conclusion section, "The pyrrolidine fragment was held well by Asp32 or Asp32′ from the native binding pose of the two compounds. Thus the interactions of other fragments with other residues in different protein regions caused significant differences." Moreover, "His66′ with its cation–pi interaction with an aniline ring of inhibitor was a perfect supporter to Asp32′. This residue's effect was absent in the other form of protonation state and caused a significant difference in simulation time. In the weak inhibitor unbinding pathways, Trp98 and Trp98′ with pi-pi interactions, due to their close position to one of the exit areas, were not good supporters for Asp32 or Asp32′, like His66′. His66′, due to far position from the bottom and center of the binding pocket, could fixed aniline fragment and decreased ligand fluctuations. The two mentioned tryptophan were closer to the important aspartic acids, and there was enough space for ligand fluctuations. For this reason, Asp36 in the active site that was close to the exit area could be a competitor with Asp32′ and was not a good interaction for keeping the ligand in the binding pocket. Similarly, attenuating effect of Trp98/Trp98′ residues in unbinding pathways of the weak inhibitor correlates with another research result" https://doi.org/10.1371/journal.pone.0257916

Comment #8:

In the main text it was mentioned that the figure 4 c-e was the first state of rep1, 2, and 3, but in figure it was mentioned as ligand RMSD. Make it correct.

Answer #8:

In Figure 4 , ligand RMSD values are for 4a and 4b.

Comment #9:

In figure 5b and c, both lines were depicted as Rep2 and Rep3-aniline-benzene. Then how does it show the difference in the distance? What are the main differences noticed in these structures?

Answer #9:

These plots aim to show and prove pi-pi intramolecular interactions in the inhibitors. So we considered the two aniline-benzene rings in compound 10 and the two nitrobenzene-benzene rings in compound 9, which were closer together to make intramolecular interactions. These plots should not compare with each other.

Comment #10:

Is the author using compound 9 as the weaker inhibitor in this study?

Answer #10:

Yes

Reviewer #3:

Comment #1:

After lots of simulations and analyses, it is still unclear that which states of protonation plays the deciding role there. It would be better if the authors are critical on this point because it is one of the most important viewpoints of this work.

Answer #1:

According to all forms of protonation, there are 16 forms. A QM/MM study will help us determine which protonation forms are closer to the natural form. Also, Aspartic proteases tend to cleave substrate by monoprotonation of aspartates in the catalytic dyad. Since HIV protease is an aspartic protease, and unfortunately, there was no typical study about HTLV-1 protease, we decided to take a pattern from HIV research. So in some research, they do not apply protonation or do not mention it, or in some research, they apply protonation state in chain A. Due to we wanted to have a complete outcome, we decided to use protonation in both chains separately.

Comment #2:

Authors in many places discussed about pi-cation and pi-anion interaction (e.g., His66 interacts with aniline in cation-pi interaction, Asp36 forming pi-anion interaction with nitrobenzene etc.). It is important to note that the energetics of such interactions are very important. (For example, pi-anion interactions have energetics in the range around 20 to 50 kj/mol whereas pi-cation interactions have energetics in the range of -10 to -20 kj/mol). The authors should check the energetics of such interactions here or at least cite relevant papers where people have reported the energetics of these particular system.

Answer #2:

We first reported the interaction energy values, so we cannot compare them or cite other works. We chose significant interactions based on the average amount of most important interaction energies of the protein-ligand complex. Our study values are negative and favorable since the nitro group is a strong electron acceptor and positively attaches to the ring. It can easily be made by pi-anion interaction with aspartic acid, which carries a negative charge.

Comment #3:

Authors should rigorously work on the figure qualities. Those are not good at publication format presently. For example, Figure 13, c, d, e. Note that, the labellings are very bold and not as per the mark of publication. Also, in 13a, the x-axis labellings are very closely spaced (10,20,30,40….80) makes it looks worse. They should scrutinize the figure qualities and font labelling throughout in-depth.

Answer #3:

The requested changes were applied to the figures.

Comment #4:

In Figure 10 caption (where they emphasize the details of compound 9 unbinding pathway), just before the last line, they mentioned "k, l, 3", What is "3" here? I did not find anything such and believe "3" should be replace by "m". Such things they should check thoroughly.

Answer #4:

The corrections were done.

Comment #5:

It is better to follow a uniform free energy bar throughout which is currently missing. For example, Figure 7, c, d, e. Their energy bars are not uniform (For c, it is 0 to 16.90 and it is different for d and e.). The authors should follow uniform energy bar for all, otherwise the comparison is not at all logical.

Answer #5:

The FEL plots were prepared by python, and the python codes made the energy bars automatically from different values obtained using "gmx sham" files in different replicas. Also, we tried uniforming plots scale, but some of the plots became too compact and became difficult to separate steps.

Comment #6:

The authors mention that "For our purpose, we had various short and long-time-scale simulations". I understand they chose based on the unbinding time scales, still did the authors performed any convergence tests or something like that, to choose time scales in this study?

Answer #6:

We did not do that. Because all full unbound occurred automatically without applying any biasing force. Various simulation times showed various unbinding processes. We corrected the text to prevent misinterpretation ( line 208).

Comment #7:

Figure 5, for the three replica's they chose three different time scales, that should be addressed properly. What happened if for Figure 5b and 5c, they did 4000 ns simulations like Figure 5a? Then, how the distance plot would look like? Since, they are doing unbiased simulations, these questions need to be answered. Also, it is visually more sound to have similar time-scales for different replica's.

Answer #7:

The time it takes for an inhibitor to leave its binding pocket depends on the pathway it is sampled. Unbinding pathways may contain multiple energy barriers, so the residence time of the inhibitor at the binding site increases, and it leaves the protein later. The inhibitor or replica whose unbinding pathway has fewer energy barriers exits the protein faster. All replicas are simulated under identical conditions using SuMD.

Comment #8:

Figure 3g, take care of the x axis label. It has been cut-off (120 ns is cut-off there). Same goes for Figure 8c. Also, The TOC (as well as TOC fonts) is not looking good as per the journal standard.

Answer #8:

The correction was done, and also revised version of the article was changed to meet PLOS ONE's style requirements.

Comment #9:

The article is somewhat looks to be less referenced. For an example, they should cite some important references regarding the interactions such as the review articles sastry et al doi: 10.1021/cr300222d and doi: 10.1021/jp900013e, and some others, for example https://doi.org/10.1021/ct200569x,https://doi.org/10.1021/acs.jpcb.7b01736and https://doi.org/10.1021/acs.jpcb.9b06343 etc.

Answer #9:

Some references were added in lines 224 and 340.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

23 May 2022

Submitted filename: Plos_review.docx

Click here for additional data file.

9 Jun 2022

Dear Prof. Sandipan Chakraborty;

I am pleased to mail you the second revision of the article titled "Comparative study of the unbinding process of some HTLV-1 protease inhibitors using Unbiased Molecular Dynamics simulations". All comments have been answered upon reviewers' recommendations, and the corrections were made in manuscript #PONE-D-22-00753R1. It is our pleasure to hear your feedback and any other suggestions relevant to the paper.

Sincerely,

Hassan Aryapour

Comment #1:

The author mentioned to my first query that “So in some research, they do not apply protonation or do not mention it, or in some research, they apply protonation state in chain A” …. It would be good if they discuss those studies briefly in the revised manuscript to get the reader a subtle grasp about the history of the study which is still somewhat missing.

Answer #1:

The requested description was added along with the reference (lines 113-117)

Comment #2:

According to my comment 5 regarding uniform free energy bar…. I think it’s needed (Even if not in the manuscript, then in the SI) because that not only establish the fidelity of the work but also required for proper comparison between the states. This is because, we know the free energy changes dramatically based on scaling and non-uniformity of the values selected.

Answer #2:

Yes, it is definitely easier to compare between states if the maximum values on the x and y scale are the same in all graphs. So once again, according to your request, we considered the maximum values of the x and y axes to be the same. Unfortunately, since in some replicas, the amount of maxima is significantly different from each other, a lot of space on the plots becomes empty and white, so the maximum amount of axes is automatically adjusted based on the maximum amount of data after running the script. However, the changes you requested were added to the supplementary section.

Comment #3:

I am somewhat surprised regarding the revision of the author of my comment 9. They did not incorporate none of the suggested studies and just incorporated two different references. Those suggested literature evidence are important and need to be incorporated.

Answer #3:

All the requested references were added (refs : 18, 38,39).

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

16 Jun 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Sandipan Chakraborty

Academic Editor

PLOS ONE

Additional Editor Comments:

Authors successfully address the reviewer concerns. However, the following issue needs to be clarified before any final decision on the manuscript.

The author of the manuscript entitled "Comparative study of the unbinding process of some HTLV-1 protease inhibitors using Unbiased Molecular Dynamics simulations" performed ligand unbinding from the protein HTLV-1. They consider one high affinity ligand and a low affinity ligand. However, the author published a previous paper on "Comparative analysis of the unbinding pathways of antiviral drug Indinavir from HIV and HTLV1 proteases by supervised molecular dynamics simulation" where they consider unbinding of another ligand Indinavir. The method is very similar and presentation is also very similar. A small discussion is there in the manuscript. But the author needs to justify their objective to work on different ligands and clearly define the novelty of the work in light of previous publication.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

27 Jun 2022

Dear Prof. Sandipan Chakraborty;

I am pleased to mail you the second revision of the article titled "Comparative study of the unbinding process of some HTLV-1 protease inhibitors using Unbiased Molecular Dynamics simulations". All comments have been answered upon reviewers' recommendations, and the corrections were made in manuscript #PONE-D-22-00753R1. It is our pleasure to hear your feedback and any other suggestions relevant to the paper.

Sincerely,

Hassan Aryapour

Comment:

Authors successfully address the reviewer concerns. However, the following issue needs to be clarified before any final decision on the manuscript.

The author of the manuscript entitled "Comparative study of the unbinding process of some HTLV-1 protease inhibitors using Unbiased Molecular Dynamics simulations" performed ligand unbinding from the protein HTLV-1. They consider one high affinity ligand and a low affinity ligand. However, the author published a previous paper on "Comparative analysis of the unbinding pathways of antiviral drug Indinavir from HIV and HTLV1 proteases by supervised molecular dynamics simulation" where they consider unbinding of another ligand Indinavir. The method is very similar and presentation is also very similar. A small discussion is there in the manuscript. But the author needs to justify their objective to work on different ligands and clearly define the novelty of the work in light of previous publication.

Answer:

We said in lines 59-69 that many research groups had designed inhibitors, but the FDA has approved none. For this reason, if we want to design successful inhibitors, we must examine existing inhibitors' negative and positive points. These will ultimately lead to valuable information for drug design.

In a previous project, our research team discovered the reasons for Indinavir's ineffectiveness against HTLV-1 protease. As we continue, we examined the unbinding pathways of two inhibitors (the most potent and one of the weakest) in complex with HTLV-1 protease. Although the study method was the same, the aim of the research was different since, in the previous project, we had one drug and two different target proteins, while in the present study, we have two different inhibitors for only one target protein.

More information was added in the Introduction section (lines 98-100), and the common results of these two projects were mentioned before in lines 456-457.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

28 Jun 2022

Comparative study of the unbinding process of some HTLV-1 protease inhibitors using Unbiased Molecular Dynamics simulations

PONE-D-22-00753R3

Dear Dr. Aryapour,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Sandipan Chakraborty

Academic Editor

PLOS ONE

4 Jul 2022

PONE-D-22-00753R3

Comparative study of the unbinding process of some HTLV-1 protease inhibitors using Unbiased Molecular Dynamics simulations

Dear Dr. Aryapour:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Sandipan Chakraborty

Academic Editor

PLOS ONE