Reconstruction of the unbinding pathways of noncovalent SARS-CoV and SARS-CoV-2 3CLpro inhibitors using unbiased molecular dynamics simulations.

The main protease (3CLpro) is one of the essential components of the SARS-CoVs viral life cycle, which makes it an interesting target for overpowering these viruses. Although many covalent and noncovalent inhibitors have been designed to inhibit this molecular target, none have gained FDA approval as a drug. Because of the high rate of COVID-19 pandemic development, in addition to laboratory research, we require in silico methods to accelerate rational drug design. The unbinding pathways of two SARS-CoV and SARS-CoV-2 3CLpro noncovalent inhibitors with the PDB IDs: 3V3M, 4MDS, 6W63, 5RF7 were explored from a comparative perspective using unbiased molecular dynamics (UMD) simulations. We uncovered common weak points for selected inhibitors that could not interact significantly with a binding pocket at specific residues by all their fragments. So water molecules entered the free binding S regions and weakened protein-inhibitor fundamental interactions gradually. N142, G143, and H163 are the essential residues, which cause key protein-ligand interactions in the binding pocket. We believe that these results will help design new potent inhibitors against SARS-CoV-2.

Introduction

Severe Acute Respiratory Syndrome (SARS) occurred in Guangdong Province of China in 2002–2003, which was caused by SARS-CoV-1, a coronavirus of 2b β-coronavirus [1]. A novel coronavirus (2019-nCoV) was identified in Wuhan, China, in late 2019 for the first time [2]. This virus, which was scientifically named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has infected many people in different parts of the world with its high prevalence power and caused the COVID-19 pandemic with symptoms including fever, cough, and fatigue [3, 4]. Sadly, now after two years, according to the World Health Organization(WHO) reports, 4,777,503 people have lost their lives all around the world (Updated on October 03, 2021) [5].

Coronaviruses are single-stranded positive-sense RNA viruses with the largest genome, approximately 30 kilobases, among all known RNA viruses. In all Coronaviruses, the genome expression is encoded by the open reading frame (ORF) 1a/b at the 5’ end of the genome [6]. Studies have shown that SARS-CoV-2 genes possess about 80% nucleotide identity and 89.10% nucleotide similarity with SARS-CoV genes [7]. Among CoVs, the viral genome of SARS-CoV-2 is about 29.8 kilobase, which encodes two polyproteins that are responsible for viral replication and transcription by proteolytic processing. This virus needs two cysteine proteases for these processes, papain-like protease (PLpro) and main proteinase (Mpro), which is also called 3C-like proteinase (3CLpro). Based on pairwise sequence alignment, 96.08% and 98.7% identity and similarity were observed between 3CLpro of SARS-CoV and SARS-CoV-2, respectively [8]. The main protease is a dimer protein with 306 residues and three domains. The domains I and II have antiparallel β barrel structures responsible for the catalytic reaction, while domain III has α-helices and regulates dimerization of the 3CLpro (Fig 1A and 1B). Since only the dimeric form of Mpro is catalytically active, intermolecular interactions between the helical domains play an essential role in activating the enzyme [9, 10].

Fig 1

The 3D ribbon structure of the 3CLpro (PDB ID 3V3M) complexed with ML188 A, The protein domains: Domain I (r8-r101), Domain II (r102-r184), Domain III (r201-r306). B, The important active site residues.

So, this protein is essential for the viral life cycle and is an interesting target for designing SARS-preventing drugs. To this end, various research groups have been working worldwide, and various approaches were used, including drug repurposing, structure-based design, and fragment-based design [11]. Based on the substrate specificity of 3CLpro, peptidomimetic inhibitors were designed as the first protease inhibitors generation [12]. These inhibitors contain Michael acceptors, aldehydes, epoxy-ketones, halo-methyl, and trifluoromethyl ketones, and they form a covalent bond with the catalytic Cys145. In the following, these inhibitors shed light on the idea for further inhibitor design, like nonpeptidic inhibitors with different micromolar ranges [13, 14]. For the SARS-CoV 3CLpro, there are some covalent and noncovalent inhibitors. At first glance, covalent warheads may seem a priority for overpowering this cysteine protease, but toxicity is one of the major challenges for the therapeutic use of these inhibitors [15]. For this purpose, we focused on reversible noncovalent inhibitors, selected four different inhibitors according to their structures and biological activity. Two of the selected compounds are SARS-CoV 3CLpro inhibitors (PDB IDs: 3V3M and 4MDS), and the other two are SARS-CoV-2 3CLpro inhibitors (PDB IDs: 6W63 and 5RF7) (Fig 2A–2D).

Fig 2

The fragmented form of the 2D structure of selected inhibitors and the.

A, The structure of ML188 in PDB ID: 3V3M and occupied positions of binding pocket by inhibitor are represented by "P" letters. B, The structure of ML300 in PDB ID: 4MDS and occupied positions of binding pocket by inhibitor are represented by "P" letters. C, The structure of inhibitor3 in PDB ID: 6W63 and occupied positions of binding pocket by inhibitor are represented by "P" letters. D, The structure of inhibitor4 in PDB ID: 5RF7 and occupied positions of binding pocket by inhibitor are represented by "P" letters.

Unfortunately, none of the designed inhibitors has been approved by the FDA. So it is pretty clear that along with experimental researches, in silico methods like unbiased molecular dynamics (UMD) are essential [16]. MD simulation has been used in many science fields since the 1950s to predict hidden information that cannot be reached through experimental research [17]. The UMD method eliminates artificial interactions between protein and inhibitor because it does not apply any biasing forces or potentials to the simulations [18]. Studying the unbinding mechanisms of inhibitors in complex with their target proteins is one of the great features of the UMD simulation method. In judging candidate drugs, the bioavailability, selectivity, metabolic properties, and binding affinity of the designed inhibitor to its target protein are important [19]. In addition to these parameters, the mean lifetime that the drug remains in the binding site is equally important. Experimental techniques can measure the time it takes for a drug to unbind from a target, but the essence of the matter is much deeper than a number [20]. On the other hand, by investigating the unbinding pathways of particular inhibitors, important information involving: protein-ligand key interactions, ligands interacting efficiently with the target can be obtained. Finally, a fully atomistic scenario will be presented based on the obtained results [21]. As a result, many research groups have examined unbinding pathways of various drugs or inhibitors over the years via MD simulation methods and prepared a solid foundation for rational drug design. [22, 23]. Among advanced MD approaches, the supervised molecular dynamics (SuMD) [24, 25] method is relatively novel. With a tabu-like supervision algorithm, it is possible to fully unbind small molecules from their molecular targets within very short times without applying any biasing force or potential. This method also produces information regarding metastable intermediate ligand-bound states, which are essential for rational drug design. In this regard, the SuMD was used to examine different cases of protein-ligand recognition mechanism, involving: the human casein kinase 2 (CK2) complexed with ellagic acid, the P1-1 isoform of glutathione S-transferase (GSTP1-1) complexed with sulfasalazine, the human peroxiredoxin 5 (PRDX5) complexed with benzene-1,2-diol, and the human serum albumin (HSA) in complex with (S)-naproxen [25].

Half-maximal inhibitory concentration (IC50) of SARS-CoV 3CLpro inhibitors are measured experimentally before, and the 3CLpro of SARS-CoV and SARS-CoV-2 are approximately the same, without any difference in the binding site [14, 26]. Therefore, based on participation in the solidarity clinical trial of COVID-19 treatments, our research team decided to compare the noncovalent SARS-CoV and SARS-CoV-2 protease inhibitors’ unbinding pathways using the SuMD.

Methods

All simulations were originated from X-ray crystallography of 3CLpro-ligand complex in Protein Data Bank (PDB IDs: 3V3M [14], 4MDS [26], 6W63 [27], 5RF7 [28]). At first, missing atoms and residues of proteins were added and fixed using UCSF Chimera software [29]. Then ligands were parameterized by ACEPYPE using default settings (the GAFF atom type and BCC partial charges) [30]. After preparation, protein-ligand complexes were constructed in GROMACS 2018 [31] using AMBER99SB force field [32] and TIP3P water model [33]. Selected holo-proteins were located in the center of triclinic boxes with a distance of 1.2 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. Using the Linear Constraint Solver (LINCS) algorithm, all Covalent bonds were constrained to maintain constant bond lengths [34]. The long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME) method [35], and the cut-off radii for Coulomb and Van der Waals (VdW) short-range interactions were set to 0.9 nm for all systems. Finally, the modified Berendsen (V-rescale) thermostat [36] and Parrinello-Rahman barostat [37] were applied for 100 and 300 ps for the equilibrations and keep the system in stable environmental conditions (310 K, 1 Bar) and got ready to begin molecular dynamic simulations with a time step of 2 fs and without applying any human or non-human biasing force or potential. In this regard, to reach complete unbinds, we performed 12 separate series of replicas (three replicas for each complex), with fixed duration times by the SuMD method with some modifications. Herein, we set the center of mass (COM) of ligands as a first spot, and the COM of His41, Cys145, His163, Asp187 in PDB ID 3V3M, His41, Met49, Cys145, His164 in PDB ID 4MDS, and His41, Cys145, His163, Met165, Gln189 in PDB IDs 6W63 and 5RF7 as second spots and, ran all simulations with a time window of 500 ps. After finishing each run, the frame with the longest distance between selected spots was selected automatically to extend the next 500 ps simulation. These processes were continued until complete unbinding 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. Figures were created using UCSF Chimera and Daniel’s XL Toolbox (v7.3.4) [38]. In addition, Matplotlib was used to create the free energy landscape plots to visualize the essential interactions [39]. The free energy landscapes plots were made based 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 color indicate that the inhibitor has been present in this area for a longer time.

Results and discussion

One of the selected compounds, ML188 (PDB ID 3V3M) with an IC50 of 4.11 μM [14], were simulated in 3 replicas, and the complete unbinding processes occurred at the times of 60, 40, and 37 ns (Fig 3A). In the first and longest replica, during the first state (Fig 3K), the ligand was enclosed within the binding pocket VdW forces for 55 ns (Fig 3G). However, all residues were not equally important; the most prominent VdW interactions formed between Met49 and Met165 residues and butylphenyl fragment of the inhibitor in its binding pose (Fig 3B). By rotation of the furan ring in inhibitor, pyridine fragment formed the third prominent VdW and amino-pi interactions with Gln189 (Fig 3C and 3G). Presumably, this rotation occurred because, unlike the butylphenyl fragment, which was well held by two methionines, other fragments did not significantly interact with the binding pocket at a specific residue. So the ligand moved from the deeper part of the binding pocket toward the exit area. These interactions became weaker due to the formation of H-bonds between the oxygen atoms in the furan ring and Ala46 and Glu47’s backbones. (Fig 3F and 3H). Since these H-bonds pulled the ligand out of the catalytic site completely (Fig 3D), the last protein-ligand interactions in the second short intermediate state (~ 10 ns) cannot be considered as essential bonds because they formed out of binding pocket and just increased simulation time (S1 Video).

Fig 3

The details of ML188 unbinding pathways in three replicas.

A, RMSD values of the ligand from binding pose to complete unbinding in three replicas. B, The interactions between particular ligand fragments and essential residues in the crystallographic binding pose of rep1. C, The new interactions between the inhibitor and binding pocket residues after rotation of the furan ring in rep1 (frame 3030 in the trajectory file). D, The last protein-inhibitor interactions before complete unbinding in rep1(frame 5804 in the trajectory file). E, The competition between C145 and butylphenyl in interaction with His41 in rep 2 and 3 (frame 27 in trajectory file of rep2). F, Hydrogen bond numbers of Ala 46, Glu47 with furan fragment in the second intermediate state of rep1. G and H, The average of most important interaction energies of the protein-ligand complex in the first and second intermediate state of rep1, respectively. I and J, The average of most important interaction energies of the protein-ligand complex in rep2 and 3, respectively. K, L and M, The free energy landscape of rep1, 2, and 3 to capture lowest energy stable states of ligand during the unbinding process (bound state (B), intermediate state (I), unbound (U)), respectively, which was calculated using "gmx sham".

In the two other replicas (rep2 and 3), the inhibitor unbound completely sooner due to the lack of H-bond formation with Ala46 and Glu47, but in comparison with the rep1, in addition to those three essential residues, His41 with VdW interactions with butylphenyl group of ligand was the fourth prominent residue (Fig 3E, 3I and 3J). The His41, due to its good position in the binding pocket, could play a vital keeping role in their single states (Fig 3L and 3M). However, because this residue had pi-sulfur interaction with Cys145, it sometimes loosed its effect on butylphenyl. So the inhibitor among the competition with Cys145 gradually pulled toward Met49 and Gln189 with the help of furan fragment rotations (S2 Video). In laboratory research on this inhibitor, Gly143, Cys145, and His163 were the most important inhibitor-protein interactions in the binding pose. Compared with our results during three replicas, Gly143 and His163 were less important than other critical residues, but Cys145 was important in the rep3 [14].

For this compound, due to the furan-free fragment, water molecules entered into a deep part of the binding pocket in the simulations. (Fig 4A and 4B). These molecules promoted unbinding by gradually weakening noncovalent interactions between the inhibitor and protein, thereby allowing a complete unbinding. (Fig 4C).

Fig 4

The details of solvation effects on ML188 unbinding mechanisms.

A, Number of water molecules in the cut-off of 3.5 Å of the binding pocket residues, rep1, 2, and 3. B, Number of water molecules in the cut-off of 5 Å of the inhibitor, rep1, 2, and 3. C, The total interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

In the following, the inhibitor, ML300 (PDB ID 4MDS), was selected as a second micromolar noncovalent inhibitor with an IC50 of 6.2 μM [26]. In the single state of the first replica in 22 ns (Fig 5A and 5L), pyrrole fragment of inhibitor was in VdW and pi-sulfur interactions with Met49, Met165, and amino-pi and VdW interactions with Gln189 (Fig 5B). Also, the benzotriazole fragment had VdW interactions with Asn142 (Fig 5H), and the O1 atom close to the benzotriazole fragment had H-bond interaction with the backbone of Glu166 (Fig 5G). Except for these two fragments, the other two were not in serious keeping interactions. So by time passing, acetylamino phenyl fragment, due to its close position to Asn142, entered in the competition with benzotriazole fragment and destroyed the effect of Asn142 on benzotriazole (Fig 5C). The most potent inhibitor interaction with the binding pocket was hydrogen bond with Glu166, so by weakening the VdW network (Fig 5F), the ligand unbound from the Glu166 side (S3 Video).

Fig 5

The details of ML300 unbinding pathways in three replicas.

A, RMSD values of the ligand from binding pose to complete unbinding in three replicas. B, The interactions between particular fragments of ligand and important residues in the crystallographic binding pose of rep1 C, The new interaction between the acetylamino phenyl fragment and Asn142 after its rotation in rep1(frame 1207 in the trajectory file). D, The interactions between particular fragments of ligand and important residues in rep2 (frame 372 in the trajectory file). E, The last interactions of ligand and protein in the second intermediate state of rep3 (frame 1172 in the trajectory file). F, The average of most important residues Lennard-Jones (LJ) energies of the protein-ligand complex in rep1, 2, and 3 as a function of time. G, Hydrogen bond numbers of Glu166 with the oxygen atom close to benzotriazole fragment in rep1 and 3, and Ser139 with acetylamino phenyl fragment in rep3. H and I, The average of most important interaction energies of the protein-ligand complex in rep1 and 2, respectively. J and K, The average of most important interaction energies of the protein-ligand complex in the first and second intermediate states of rep3. L, M and N, The free energy landscape of rep1, 2, and 3 during the unbinding process (bound state (B), intermediate state (I), unbound (U)), respectively.

In the other pathway (rep2), by the time of 24 ns, pyrrole fragment had VdW and pi-sulfur interactions with Met49 and Met165, and in contrast with rep1, Gln189 did not have a significant effect (Fig 5J) in the single state of this replica (Fig 5M). Also, the acetylamino phenyl group was not free to break down the interaction of the benzotriazole and Asn142 because it had VdW and cation-pi interactions with His41 (Fig 5D). Furthermore, the most important difference between rep1 and rep2 was the absence of a critical hydrogen bond at Glu166, which caused the ligand to exit the protein from its other side (S4 Video).

Ultimately, in the last and also in the shortest replica with the time of 18 ns, due to lack of significant catalytic VdW forces (Fig 5F and 5J), Glu166 by forming hydrogen bond (Fig 5G), caused the ligand to be pulled out of the active site in the first state. In the second intermediate state (Fig 5N), Tyr126 and Ser139 residues by VdW (Fig 5K) and hydrogen bonding interactions with acetylamino phenyl kept the inhibitor in the protein exposure, out of the binding pocket for 14ns, respectively (S5 Video). In the results of experimental research, a long list of important residues is reported. Between these reported important residues, Met49, Glu166, and Gln189 correlate with our results [26].

In contrast with ML188, more water molecules in the native binding pose of this inhibitor caused essential interactions to become water-mediated from the first moment of simulation (Fig 6A and 6B). Also, by time passing, due to benzotriazole free fragment’, more space was created for water molecules insertion into the binding pocket, and ligand got unbound more rapidly (Fig 6C).

Fig 6

The details of solvation effects on ML300 unbinding mechanisms.

A, Number of water molecules in the cut-off of 3.5 Å of the binding pocket residues, rep1, 2, and 3. B, Number of water molecules in the cut-off of 5 Å of the inhibitor, rep1, 2, and 3. C, The total interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

We proceeded to simulate the unbinding mechanism of two additional compounds to confirm and complete the information obtained and perhaps even identify new key factors. The potent compound in PDB ID 6W63, with a broad spectrum of anti-viral activities, was chosen to achieve this goal. Furthermore, since this compound has some structural similarities to ML188, we were interested in understanding what was happening between this inhibitor and the protein while it was unbinding. In the first and second replicas with 55 ns and 61 ns, respectively (Fig 7A), the pyridine ring of the compound had cation–pi interaction with His163 in the deep part of the binding pocket (Fig 7B), and also the backbone of Gly143 had H-bond interaction with the oxygen atom of imidazole fragment in the shallow part of the pocket (Fig 7F). These potent interactions and VdW interaction between Met165 and butylphenyl fragment were the most important protein-inhibitor interactions in the native binding pose (Fig 7H and 7I). Later, Asn142, by forming a H-bond with the oxygen atom of the imidazole fragment (Fig 7C and 7G), caused this fragment to be kept by two consecutive residues. This good binding pose was continued until 30ns, insofar as, Asn142 switched its H-bond to an imidazole ring (Fig 7D), and while Met 165 was still in interaction with the butylphenyl fragment, His163 lost its effect. So the ligand gradually moved from the deep part of the binding pocket to the surface area. Then, after time passing, by rotation of butylphenyl fragment, cyclohexylacetamide fragment formed H-bond with Glu166 (Fig 7E), and finally, the inhibitor unbound from Glu166 side.

Fig 7

The details of the inhibitor3 in PDB ID 6W63 unbinding pathways in three replicas.

A, RMSD values of the ligand from binding pose to complete unbinding in three replicas. B, The interactions between particular fragments of ligand and important residues in the crystallographic binding pose of rep1, 2, and 3. C, The interaction between the imidazole fragment and Asn142 in rep1, 2, and 3 (frame 607 in the trajectory file of rep1). D, The new interaction between Asn142 and imidazole ring in rep1, 2, and 3 (frame 2843 in the trajectory file of rep1). E, The new interaction between Glu166 and cyclohexylacetamide fragment, after butylphenyl fragment rotation in rep1, 2, and 3 (frame 4291 in the trajectory file of rep1). F, Hydrogen bond numbers of Gly143 and O1 and also Hydrogen bond numbers of Glu166 with cyclohexylacetamide in all replicas G, Hydrogen bond numbers of Asn142 and O1, and also imidazole ring in all replicas. H, I, and J, The average of most important interaction energies of the protein-ligand complex in rep1, 2, and 3, respectively. K, L, and M The free energy landscape of rep1, 2, and 3 during the unbinding process (bound state (B), intermediate state (I), unbound (U)), respectively.

In the final replica, due to the lack of a continuous binding pose key hydrogen bonds (Fig 7F and 7K), the ligand left the protein after 30 ns (Fig 7A). With more details, the lifetime of H-bond interaction between Gly143 and butylphenyl fragment was too short (95 ps), so this fragment could not be fixed and exposed to Met165. Like the last two replicas, butylphenyl fragment rotation and complete unbinding were observed sooner (Fig 7J). This inhibitor had only one state (Fig 7K–7M), and after weakening the binding pose interactions, the ligand did not trap in a serious state (S6 Video).

At the beginning of the unbinding process of the inhibitor3 in PDB ID 6W63, the crystalline water molecules are present at the native binding pose. In the following, more water molecules came into the binding site and broke important interactions (Fig 8A and 8B). Ultimately the ligand is unbound with complete solvation (Fig 8C).

Fig 8

The details of solvation effects on inhibitor3 in PDB ID 6W63 unbinding mechanisms.

A, Number of water molecules in the cut-off of 3.5 Å of the binding pocket residues, rep1, 2, and 3. B, Number of water molecules in the cut-off of 5 Å of the inhibitor, rep1, 2, and 3. C, The total interactions energies of protein-inhibitor complexes in rep1, 2, and 3.

Finally, the inhibitor in PDB ID 5RF7 was the last candidate because of its different structure compared with ML188. In this regard, all its replicas with times of 12, 16, and 17 ns, consider as a rapid unbinding pathway (Fig 9A). So in the binding pose of this pathway, while the 1,4-dimethylpiperazine fragment had only ‌pi-sulfur interaction with Met165, the pyrrolo[2,3-b] pyridine fragment of inhibitor was in VdW and amino-pi interactions with Asn142, VdW interaction with Glu166, and had polar–pi interaction with Ser144 and also had Cation–pi interaction with His163 (Fig 9E–9G). Pyrrolo [2,3-b] pyridine, unlike 1,4-dimethylpiperazine fragment, was well kept by various interactions with both superficial and deep residues (Fig 9B). In the following, by rotation of methylpiperazine fragment toward Glu166, His163 loosed its strategic effect to let the ligand enter to the second intermediate state of the two longer replicas (rep1 and rep3) and be tapped in anion-pi interaction of Glu166 with pyrrolo[2,3-b] pyridine fragment (Fig 9C) as the last protein-inhibitor interaction (S7 Video). The second intermediate state of rep1 and 3 was unstable, as was the first state of all replicas (Fig 9H–9J). Even though all these inhibitor fragments were in serious interactions and there were no water molecules between these fragments and specific residues in the binding pose, there was enough space for water molecules to enter (Fig 9K and 9L). So this tiny ligand, by water-mediated interactions, loosed all its important interactions and unbound in a quick time (Fig 9D).

Fig 9

The details of the inhibitor4 in PDB ID 5RF7 unbinding pathways in three replicas.

A, RMSD values of the ligand from binding pose to complete unbinding in three replicas. B, The interactions between particular fragments of ligand and important residues in the crystallographic binding pose of rep1, 2 and 3. C, The new interaction between the pyrrolo[2,3-b] pyridine fragment and Glu166 in rep1, 2, and 3 (frame 548 in the trajectory file of rep1). D, The total interaction energies of protein-inhibitor complexes in rep1, 2, and 3. E, F, and G The average of most important interaction energies of the protein-ligand complex in rep1, 2, and 3, respectively. H, I, and J, The free energy landscape of rep1, 2, and 3 during the unbinding process (bound state (B), intermediate state (I), unbound (U)), respectively. K, Number of water molecules in the cut-off of 3.5 Å of the binding pocket residues, rep1, 2, and 3. L, Number of water molecules in the cut-off of 5 Å of the inhibitor, rep1, 2, and 3.

Conclusion

By putting together the atomic details of the unbinding pathways of selected inhibitors except inhibitor 4 in PDB ID 5RF7, one significant common weakness point was observed at various times in all replicas. There were no serious interactions between all fragments of inhibitors at specific residues in the binding pocket, so this factor was sufficient to weaken critical interactions. Almost free fragments could compete with other fragments for interactions with key residues. Even if they did not engage in competitive interactions, they could still change the inhibitors’ positions and move closer to the exit path by rotating them. So the inhibitor in PDB ID 6W63 made serious noncovalent interactions with N142, G143, and H163 in the binding pocket, but due to its free 1-(1H-imidazol-4-yl) ethanone fragment could not stay longer in the binding pose.

In another study, hydrogen bonding interaction between inhibitor in PDB ID 3V3M and His163 was considered an essential interaction, but in our three replicas, the role of His163 in unbinding was not critical [14]. Based on the other research for the inhibitor in PDB ID 4MDS, Met49 and Gln189 were important residues correlated with our important residue list in keeping the inhibitor at the binding site [26].

Furthermore, water molecules played a functional role in all unbinding mechanisms by interfering and breaking important protein-inhibitor interactions. As time progressed, all inhibitors could not interact with all S regions of the binding pocket, and there was enough space for more water molecules to be inserted from outside. So the inhibitor 4 in PDB ID 5RF7, which did not have a free fragment, was too tiny and, in comparison with other inhibitors, occupied less space of binding pocket. So there was more space for water molecules to enter and caused less time to unbind due to the solvation effect. On the other hand, when the inhibitor occupies all space of the binding pocket, water molecules cannot penetrate under the ligand, and important inhibitor-protein interactions do not become water-mediated. In conclusion, the next series of noncovalent inhibitors should be designed to occupy all S regions of the binding pocket to make maximum noncovalent interactions. This information is valuable for designing a new generation of inhibitors against this molecular target by fixing the weaknesses mentioned.

Video1, unbinding pathway of ML188 in rep1.

(MP4)

Click here for additional data file.

Video2, unbinding pathway of ML188 in rep2 and 3.

(MP4)

Click here for additional data file.

Video3, unbinding pathway of ML300 in rep1.

(MP4)

Click here for additional data file.

Video4, unbinding pathway of ML300 in rep2.

(MP4)

Click here for additional data file.

Video5, unbinding pathway of ML300 in rep3.

(MP4)

Click here for additional data file.

Video6, unbinding pathway of inhibitor 3 in PDB ID 6W63 in rep1-3.

(MP4)

Click here for additional data file.

Video7, unbinding pathway of inhibitor 4 in PDB ID 5RF7 in rep1-3.

(MP4)

Click here for additional data file.

9 Nov 2021

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,

Chandrabose Selvaraj, Ph.D.

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. Please update your submission to use the PLOS LaTeX template. The template and more information on our requirements for LaTeX submissions can be found at http://journals.plos.org/plosone/s/latex.

3. We note that the grant information you provided in the ‘Funding Information’ and ‘Financial Disclosure’ sections do not match.

When you resubmit, please ensure that you provide the correct grant numbers for the awards you received for your study in the ‘Funding Information’ section.

4. Thank you for stating the following financial disclosure:

"YES

grant number 99-213-1 from Golestan University, Gorgan, Iran."

Please state what role the funders took in the study.  If the funders had no role, please state: "The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript."

If this statement is not correct you must amend it as needed.

Please include this amended Role of Funder statement in your cover letter; we will change the online submission form on your behalf.

Reviewers' comments:

Reviewer's Responses to Questions

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: In the article entitled “Reconstruction of the unbinding pathways of noncovalent SARS-COV and SARSCOV-2 3CLpro inhibitors using Unbiased Molecular Dynamics” the authors performed an in silico analysis of the unbinding process for two SARS-CoV and two SARS-CoV-2 3CLpro inhibitors. While the idea of the research seems good and very useful and the research seems to be performed correctly, the discussion is lacking; there are several issues that require attention and need to be corrected and more thoroughly explained before the article is accepted:

1. Capitalization of abbreviations “SARS-CoV” and “SARS-CoV-2” needs to be corrected in the title, as well as in the text. Additionally, capitalization of other words, such as “viruses” in the abstract needs to be corrected.

2. Also, there seems to be an issue with font type and size which is inconsistent.

3. Throughout the manuscript, English should be corrected e.g. the sentence (lines 35-36): “Coronaviruses (CoVs), which have the largest genome among RNA viruses, are enveloped by single-stranded RNA viruses.” does not make any sense.

4. The sentence (lines 71-72): “On the other hand, investigation of the unbinding pathways of particular inhibitors, protein-ligand key interactions, and ligand's effectiveness or weakness can be obtained.” needs to be rephrased, it does not make a lot of sense.

5. The sentence (line 73): “Ultimately, the obtained results will be presented as a fully atomistic scenario [17].” needs to be rephrased (and font changed).

6. The sentence (lines 86-90): “Since, half-maximal inhibitory concentration (IC50) of SARS-COV 3CLpro inhibitors are measured experimentally before, and the 3CLpro of SARS-COV and SARS-CoV-2 are approximately the same, our research team for participating in the global solidarity trial decided to examine some noncovalent SARS-CoV and SARS-CoV-2 protease inhibitors' unbinding pathways by the SuMD, Comparatively.” is very confusing and needs to be rephrased. Also, what does “approximately the same” even mean? Additional clarification is needed and does this influence the binding site?

7. There is no discussion and comparison of the obtained results with results of other papers, have other researchers performed similar studies, are the conclusions the same, is there any difference in results or a new insight?

8. What are the differences between the SARS-CoV and SARS-CoV-2 3CLpro binding sites? How does this reflect on binding of the studied inhibitors and their selectivity? If there are not any differences, it should also be said.

9. In figures, such as Figure 3 l, m and n, should the word “state (S)” be “bound state”?

10. Figures in general are of low quality, I don’t know if this is a consequence of converting to the PDF format or not.

11. Can the authors explain the “free energy” term in e.g. Figure 2 m, n, o, is it the ∆GBIND or something else?

12. Can the authors also explain how come in some cases (e.g. Figure 3.c) the total interaction energy goes as low as -350 kJ/mol, what does this interaction energy represent?

13. In line 255, what does the expression “the almost deep part” mean?

14. In the Conclusion, English is especially poor: the sentence (lines 336-338): “All fragments of inhibitors did not interact seriously at specific residues in the binding pocket, so this factor was sufficient to weaken critical interactions, resulting in rapid unbinding.” is very confusing and has no meaning.

15. The same goes for the sentence after that (lines 338-341): “These almost free fragments could compete with other fragments for interactions with essential residues. Even if they did not enter into competitive interactions, they could change inhibitors' position and brought the inhibitor closer to the exit path by their rotations.” which needs to be rephrased.

16. The sentence (lines 353-355): “In conclusion, the next series of noncovalent inhibitors should be designed so that they can occupy all S regions of the binding pocket to make maximum noncovalent interactions.” Does not say anything new, this is common knowledge – the inhibitors should fill out the binding site as much as possible. Are there any concrete suggestions the authors have for development of future 3CLpro inhibitors, which interactions are good and should stay, which interactions are bad and should be corrected etc.?

Reviewer #2: The article "Reconstruction of the unbinding pathways of noncovalent SARS-COV and SARS-COV-2 3CLpro inhibitors using Unbiased Molecular Dynamics" is an interesting paper in which pure classical Molecular Dynamics Simulations have been employed to investigate the unbinding process pathways of two noncovalent inhibitors of SARS-COV and SARS-COV-2. Authors have used the PDB IDs: 3V3M, 4MDS, 6W63 and 5RF7 and they have explored the inhibitors comparing the protein-inhibitor interaction by means of unbiased molecular dynamics (UMD) simulations. Authors have found that when the inhibitor occupies all space of the binding pocket, water molecules cannot penetrate under the ligand, and important inhibitor-protein interactions do not become water-mediated. In other words, the inhibitor must occupy all S region to develop a proper binding, otherwise, water molecules can occupy this space and important interactions can get lost.

Comments:

1) The methodology of this work is correct, the analysis is well supported, and the results are partially interesting.

2) The SuMD MD simulation is interesting to unbind the inhibitor from the active site because the center of masses distance is increased in each window. Perhaps, for pure classical Molecular Dynamic Simulations, just 500 ps could be not enough to equilibrate the system. Have the authors tried to increase the simulation time by window to avoid possible equilibration problems?

3) The structure of the paper and the results are very clear. There are four systems: two inhibitors and two proteins (3CL of SAR-CoV and SARS-CoV-2). For each system, three replicas unbinding the inhibitor (from 25 to 65 ns). For the results, same figures and same analysis to compare the results properly. My major criticism about this paper is that the main goals from the results: on one hand, when the inhibitor is unbinding, water molecules are coming in the active site, so, when the inhibitor occupy the S region, the binding should be better because less water molecules are affecting the inhibitor-enzyme interaction, so the main idea is that “water molecules played a functional role in all unbinding mechanisms by interfering and breaking important protein-inhibitor interactions”, which is something very well known in from many other works On the other hand, the inhibitor interacts such as independent fragments with different parts of the proteins. I think these findings are not very novel, however, some results are very interesting such as Figures 2, 4 and 6, parts m, n and o. These maps should be explained by the authors more carefully.

4) The bibliography is insufficient. There are many QM/MM Molecular Dynamic simulations and classical simulations of many inhibitors of the 3CL enzyme of SAR-CoV and SARS-CoV-2 which must be included by the authors.

5) As minor points:

1) The structure of the inhibitors should be provided in the introduction section. In addition, the structure of the active site of the complex protein-inhibitor highlighting the critical distances should be showed as well.

2) Page 9, line 178 “j and k” looks like a typo mistake

3) Footnote of Figure 3. Cut-offs are 0.35 and 0.5 angstroms? I think there is a mistake is those data.

Reviewer #3: The manuscript ID PONE-D-21-32120 entitled “Reconstruction of the unbinding pathways of noncovalent SARS-COV and SARS-COV-2 3CLpro inhibitors using Unbiased Molecular Dynamics” is quite a good study. The main protease (3CLpro) is one of the essential components of the SARS-COVs viral life cycle, which makes it an interesting target for overpowering Viruses. Although many covalent and noncovalent inhibitors have been designed to inhibit this molecular target, none have gained FDA approval as a drug. Because of the high rate of COVID-19 pandemic development, in addition to laboratory research, we require in silico methods to accelerate rational drug design. The unbinding pathways of two SARS-COV and SARS-COV-2 noncovalent inhibitors with the PDB IDs: 3V3M, 4MDS, 6W63, 5RF7 were explored from a comparative perspective using the unbiased molecular dynamics (UMD) simulations. We uncovered common weak points for selected inhibitors that could not have significant interactions with a binding pocket at specific residues by all their fragments. So water molecules entered the free binding S regions and weakened protein-inhibitor key interactions gradually. Finally, the authors believed that these results will help to design new potent inhibitors against SARS-CoV-2.

I appreciate the authors for their great effort of enclosing all the supplementary materials. The videos are really nice.

However, the following queries should be addressed before submitting a revision

1) In the abstract, the author provided more general info instead of the results of the study. In addition, the conclusion and methods are not appropriate---make it clear

2) In the introduction section, the author might start with SARS first, then SARS-2 to follow the uniformity of the title of the study.

3) Page 8 – lines 46-49: “The domains I and II have 46 antiparallel β barrel structures responsible for the catalytic reaction, while domain III has α-helices (Fig 1a, b). Since the dimeric form of Mpro is catalytically active, intermolecular interactions between the helical domains play an essential role in activating the enzyme [8, 9]. ------ It is not clear. What is the role of domain III?

4) What is unbiased molecular dynamics (UMD) and its advantages? Brief up 1 or 2 lines.

5) This line seems to be meaningless – “MD simulation has been a standard method of predicting hidden information in many science fields since the 1950s”----what is meant by hidden info by MD simulations?

6) What is the correlation of the above lines and “in judging candidate drugs, the bioavailability, selectivity, metabolic properties, and binding affinity of the designed inhibitor to its target protein are all important? [16].”

7) It is not an appropriate word “ligand's effectiveness or weakness can be obtained”. It may be ligands interacting efficiently with the target

8) These lines are a bit confusing or meaningless “Since, half-maximal inhibitory concentration (IC50) of SARS-COV 3CLpro inhibitors are measured experimentally before, and the 3CLpro of SARS-COV and SARS CoV-2 are approximately the same, our research team for participating in the global solidarity trial decided to examine some noncovalent SARS-CoV and SARS-CoV-2 protease inhibitors' unbinding pathways by the SuMD, Comparatively.” From MD authors defining IC50 value?

9) This is unclear: we performed 12 simulations, three replicas for each complex, by the SuMD method with some modifications---what are 12 simulations? Whether author’s starts from 1ps to 500ps or 1ps to 1500 ps defined as replicas? Because the following line authors mentioned, “------to extend the next 500 ps simulation”.

10) Reference missing for UCSF Chimera

11) Author may improve the discussion because in the results and discussion section authors described only results.

12) How do the authors create a free energy plot? Which package?

13) In all the figures, authors may reduce the thickness of the lines to get better visualization---figures are not in good resolution

14) Conclusion should be improved.

15) The language of the manuscript should meet the journal’s adequate standard

[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.

11 Dec 2021

Dear Prof. Chandrabose Selvaraj;

I am pleased to mail you the revised article titled "Reconstruction of the unbinding pathways of noncovalent SARS-CoV and SARS-CoV-2 3CLpro inhibitors using Unbiased Molecular Dynamics". All comments have been answered upon reviewers' recommendations, and the corrections were made in manuscript #PONE-D-21-32120. It is our pleasure to hear your feedback and any other suggestions relevant to the paper.

Sincerely,

Hassan Aryapour

Reviewer #1:

Comment #1:

Capitalization of abbreviations "SARS-CoV" and "SARS-CoV-2" needs to be corrected in the title, as well as in the text. Additionally, capitalization of other words, such as "viruses" in the abstract needs to be corrected.

Answer #1:

The requested items were corrected.

Comment #2:

Also, there seems to be an issue with font type and size which is inconsistent.

Answer #2:

The entire text was reviewed and corrected.

Comment #3:

Throughout the manuscript, English should be corrected e.g. the sentence (lines 35-36): "Coronaviruses (CoVs), which have the largest genome among RNA viruses, are enveloped by single-stranded RNA viruses." does not make any sense.

Answer #3:

A new version of this sentence is highlighted in (lines 40-42).

Comment #4:

The sentence (lines 71-72): "On the other hand, investigation of the unbinding pathways of particular inhibitors, protein-ligand key interactions, and ligand's effectiveness or weakness can be obtained." needs to be rephrased, it does not make a lot of sense.

Answer #4:

A new version of this sentence is highlighted in (lines 85-87).

Comment #5:

The sentence (line 73): "Ultimately, the obtained results will be presented as a fully atomistic scenario [17]." needs to be rephrased (and font changed).

Answer #5:

A new version of this sentence is highlighted in (lines 88-89).

Comment #6:

The sentence (lines 86-90): "Since, half-maximal inhibitory concentration (IC50) of SARS-COV 3CLpro inhibitors are measured experimentally before, and the 3CLpro of SARS-COV and SARS-CoV-2 are approximately the same, our research team for participating in the global solidarity trial decided to examine some noncovalent SARS-CoV and SARS-CoV-2 protease inhibitors' unbinding pathways by the SuMD, Comparatively." is very confusing and needs to be rephrased. Also, what does "approximately the same" even mean? Additional clarification is needed and does this influence the binding site?

Answer #6:

This paragraph was rewritten for better understanding. We used "approximately the same", based on multiple sequence alignment, which reported 96.08% identity and 98.7% similarity between 3CLpro of SARS-CoV and SARS-CoV-2. There is no difference for the 3CLpro binding sites of these two viruses, so we added this information to the new version of this paragraph.

Comment #7:

There is no discussion and comparison of the obtained results with results of other papers, have other researchers performed similar studies, are the conclusions the same, is there any difference in results or a new insight?

Answer #7:

PDB IDs 3V3M and 4MDS have other research, and we added a comparison of their results with ours in the result and conclusion sections in the (lines 187-190 & 240-242 & 371-375) but for the two other PDB IDs (6W63 and 5RF7) there are no articles.

Comment #8:

What are the differences between the SARS-CoV and SARS-CoV-2 3CLpro binding sites? How does this reflect on binding of the studied inhibitors and their selectivity? If there are not any differences, it should also be said.

Answer #8:

The comment was answered alongside Comment #6.

Comment #9:

In figures, such as Figure 3 l, m and n, should the word "state (S)" be "bound state"?

Answer #9:

In all figures, state words were replaced with bound states (B).

Comment #10:

Figures in general are of low quality, I don't know if this is a consequence of converting to the PDF format or not.

Answer #10:

We fixed the problem of low-quality figures caused by PDF conversion.

Comment #11:

Can the authors explain the "free energy" term in e.g. Figure 2 m, n, o, is it the ∆GBIND or something else?

Answer #11:

The free energy landscape plots were obtained using "gmx sham" as we mentioned in the figures. Actually, "gmx sham" leads to plot Gibbs's free energy landscapes.

Comment #12:

Can the authors also explain how come in some cases (e.g. Figure 3.c) the total interaction energy goes as low as -350 kJ/mol, what does this interaction energy represent?

Answer #12:

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. Also, -350 kJ / mol values in some plots and -300 kJ / mol values in others indicate the value of bound state interaction energies. Since these interactions measure in vacuum conditions, these values cannot prove that the inhibitor with a lower interaction energy value in a bound state is more potent. These plots should be examined along with the average of the most important interaction energies for better understanding.

Comment #13:

In line 255, what does the expression "the almost deep part" mean?

Answer #13:

The word "almost" is deleted.

Comment #14:

In the Conclusion, English is especially poor: the sentence (lines 336-338): "All fragments of inhibitors did not interact seriously at specific residues in the binding pocket, so this factor was sufficient to weaken critical interactions, resulting in rapid unbinding." is very confusing and has no meaning.

Answer #14:

A new version of this sentence is highlighted in (lines 358-361).

Comment #15:

The same goes for the sentence after that (lines 338-341): "These almost free fragments could compete with other fragments for interactions with essential residues. Even if they did not enter into competitive interactions, they could change inhibitors' position and brought the inhibitor closer to the exit path by their rotations." which needs to be rephrased.

Answer #15:

A new version of this sentence is highlighted in (lines 361-365).

Comment #16:

The sentence (lines 353-355): "In conclusion, the next series of noncovalent inhibitors should be designed so that they can occupy all S regions of the binding pocket to make maximum noncovalent interactions." Does not say anything new, this is common knowledge – the inhibitors should fill out the binding site as much as possible. Are there any concrete suggestions the authors have for development of future 3CLpro inhibitors, which interactions are good and should stay, which interactions are bad and should be corrected etc.?

Answer #16:

We explained everything, even common knowledge, because we wanted to report all atomistic details. Moreover, we explained all favorable and unfavorable interactions in the results and discussion sections for all inhibitors' unbinding processes. The most potent interactions are also listed in the conclusion section (lines 367-370).

Reviewer #2:

Comment #1:

The methodology of this work is correct, the analysis is well supported, and the results are partially interesting.

Comment #2:

The SuMD MD simulation is interesting to unbind the inhibitor from the active site because the center of masses distance is increased in each window. Perhaps, for pure classical Molecular Dynamic Simulations, just 500 ps could be not enough to equilibrate the system. Have the authors tried to increase the simulation time by window to avoid possible equilibration problems?

Answer #2:

Equilibration was applied by performing 100 ps of NVT and 300 of NPT. The 500 ps simulation is not for equilibration. The method section has been improved for better understanding.

Comment #3:

The structure of the paper and the results are very clear. There are four systems: two inhibitors and two proteins (3CL of SAR-CoV and SARS-CoV-2). For each system, three replicas unbinding the inhibitor (from 25 to 65 ns). For the results, same figures and same analysis to compare the results properly. My major criticism about this paper is that the main goals from the results: on one hand, when the inhibitor is unbinding, water molecules are coming in the active site, so, when the inhibitor occupy the S region, the binding should be better because less water molecules are affecting the inhibitor-enzyme interaction, so the main idea is that "water molecules played a functional role in all unbinding mechanisms by interfering and breaking important protein-inhibitor interactions", which is something very well known in from many other works On the other hand, the inhibitor interacts such as independent fragments with different parts of the proteins. I think these findings are not very novel, however, some results are very interesting such as Figures 2, 4 and 6, parts m, n and o. These maps should be explained by the authors more carefully.

Answer #3:

For the "On the other hand, the inhibitor interacts such as independent fragments with different parts of the proteins. I think these findings are not very novel" we explained in the answer of Comment #16 of reviewe1. Also, we added a more detailed explanation of the Free energy landscape plots in the method section.

Comment #4:

The bibliography is insufficient. There are many QM/MM Molecular Dynamic simulations and classical simulations of many inhibitors of the 3CL enzyme of SAR-CoV and SARS-CoV-2 which must be included by the authors.

Answer #4:

The QM/MM Molecular Dynamic simulations are a hybrid method that is entirely unrelated to our work. Our study is MM Molecular Dynamic simulations, which there is no related project to mention.

Comment #5-1:

The structure of the inhibitors should be provided in the introduction section. In addition, the structure of the active site of the complex protein-inhibitor highlighting the critical distances should be showed as well.

Answer #5-1:

All fragments of each inhibitor are existence in the figures of the result section. If we change the location of these pictures, the readership of this article may get board to switch to the introduction section and again back to the result section. The active site domain of the protein is explained in Figure1, and a binding pocket of the complex protein-inhibitor is picturizing separately in Figures 2, 4, 6, and 8.

Comment #5-2:

Page 9, line 178 "j and k" looks like a typo mistake

Answer #5-2:

This line was rechecked.

Comment #5-3:

Footnote of Figure 3. Cut-offs are 0.35 and 0.5 angstroms? I think there is a mistake is those data.

Answer #5-3:

True values of Cut-offs were replaced in figures' captions.

Reviewer #3:

Comment #1:

In the abstract, the author provided more general info instead of the results of the study. In addition, the conclusion and methods are not appropriate---make it clear.

Answer #1:

More information about the conclusion is added in the abstract section. Also, the Conclusion and method sections were improved.

Comment #2:

In the introduction section, the author might start with SARS first, then SARS-2 to follow the uniformity of the title of the study.

Answer #2:

The requested information was added in the introduction section (lines 31-32)

Comment #3:

Page 8 – lines 46-49: "The domains I and II have 46 antiparallel β barrel structures responsible for the catalytic reaction, while domain III has α-helices (Fig 1a, b). Since the dimeric form of Mpro is catalytically active, intermolecular interactions between the helical domains play an essential role in activating the enzyme [8, 9]. ------ It is not clear. What is the role of domain III?

Answer #3:

The requested information was added in the introduction section (line 53)

Comment #4:

What is unbiased molecular dynamics (UMD) and its advantages? Brief up 1 or 2 lines.

Answer #4:

The requested information was added in the introduction section (lines 75-77)

Comment #5:

This line seems to be meaningless – "MD simulation has been a standard method of predicting hidden information in many science fields since the 1950s"----what is meant by hidden info by MD simulations?

Answer #5:

Since the experimental research gives information based on x-ray crystallographic binding pose, results are statics and limited, but MD simulations are dynamic and give information about what is happening during unbinding pathways. Therefore, this information cannot be obtained through experimental research and is regarded as hidden.

Comment #6:

What is the correlation of the above lines and "in judging candidate drugs, the bioavailability, selectivity, metabolic properties, and binding affinity of the designed inhibitor to its target protein are all important? [16].”

Answer #6:

For better understanding, we have rewritten this sentence (lines 81-84).

Comment #7:

It is not an appropriate word "ligand's effectiveness or weakness can be obtained". It may be ligands interacting efficiently with the target

Answer #7:

The "ligands interacting efficiently with the target" was replaced with "ligand's effectiveness or weakness can be obtained" in (lines 86-87).

Comment #8:

These lines are a bit confusing or meaningless "Since, half-maximal inhibitory concentration (IC50) of SARS-COV 3CLpro inhibitors are measured experimentally before, and the 3CLpro of SARS-COV and SARS CoV-2 are approximately the same, our research team for participating in the global solidarity trial decided to examine some noncovalent SARS-CoV and SARS-CoV-2 protease inhibitors' unbinding pathways by the SuMD, Comparatively." From MD authors defining IC50 value?

Answer #8:

Experimental research provided the IC50 value. The reference is added in (line 104).

Comment #9:

This is unclear: we performed 12 simulations, three replicas for each complex, by the SuMD method with some modifications---what are 12 simulations? Whether author's starts from 1ps to 500ps or 1ps to 1500 ps defined as replicas? Because the following line authors mentioned, "------to extend the next 500 ps simulation".

Answer #9:

We improved the method section. Actually, we ran all simulations with a time window of 500 ps. We mentioned in the method section, "After finishing each run, the frame with the longest distance between selected spots 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", So each replica is a separate simulation from binding pose to complete unbinding.

Comment #10:

Reference missing for UCSF Chimera

Answer #10:

The reference is available in (line 113).

Comment #11:

Author may improve the discussion because in the results and discussion section authors described only results.

Answer #11:

The conclusion section changed to improve.

Comment #12:

How do the authors create a free energy plot? Which package?

Answer #12:

The requested information was added in the method section (lines 144-150).

Comment #13:

In all the figures, authors may reduce the thickness of the lines to get better visualization---figures are not in good resolution.

Answer #13:

All figures were rebulided again, but for more precise images, we reduced the thickness of the lines, resulting in plots that faded

Comment #14:

Conclusion should be improved.

Answer #14:

The conclusion section changed to improve.

Comment #15:

The language of the manuscript should meet the journal's adequate standard.

Answer #15:

The requested work is done.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

21 Dec 2021

9 Jan 2022

Dear Prof. Chandrabose Selvaraj;

I am pleased to mail you the revised article titled "Reconstruction of the unbinding pathways of noncovalent SARS-CoV and SARS-CoV-2 3CLpro inhibitors using Unbiased Molecular Dynamics". All comments have been answered upon reviewers' recommendations, and the corrections were made in manuscript #PONE-D-21-32120R1. It is our pleasure to hear your feedback and any other suggestions relevant to the paper.

Sincerely,

Hassan Aryapour

ACADEMIC EDITOR:

You may want to upload a letter of editing from a professional agency or native-English speaker in addition to the reviewer's comments because reviewers have expressed a concern about the language part.

Answer:

The entire document has been reviewed by someone who has a better understanding of English and checked for existing problems in the text

Journal Requirements:

Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript. If you need to cite a retracted article, indicate the article's retracted status in the References list and also include a citation and full reference for the retraction notice.

Answer:

The references list was reviewed once more carefully, and no reference was changed.

Reviewer #1:

The authors have significantly improved the manuscript and added further discussion of their results. There are still significant language issues that need to be corrected throughout the manuscript (with e.g. lines 39-40 "Coronaviruses (CoVs), an RNA virus with the largest genome, are enveloped by a single-stranded RNA virus. "– coronaviruses are enveloped RNA viruses, they are not enveloped BY a single stranded RNA virus; and lines 99-101: "So our research team, for participating in the solidarity clinical trial for COVID-19 treatments, decided to examine some noncovalent SARS-CoV and SARS-CoV-2 protease inhibitors' unbinding pathways by the SuMD, Comparatively. "which is still grammatically very incorrect). Additionally, my comment about the lack of discussion and comparison with results from other researchers is still standing. The authors have improved their discussion, but they have obtained a lot of new information which could have been more thoroughly discussed and compared, it would further increase the significance and value of their work. Nonetheless, even in the present state, I believe that this research should be accepted for publishing.

Answer #2:

The lines 39 and 40 and also lines 99-101 were replaced with the corrected version in (lines 39-40 and 111-114).

For the comparison, our research method is different from the mentioned experimental research. They examined important protein-inhibitor interactions only in a binding pose, so we could compare only their ultimate result with a little part of our work.

Reviewer #2: 1) Answer #5-1:

Reviewer #2: 1) Answer #5-1:

All fragments of each inhibitor are existence in the figures of the result section. If we change the location of these pictures, the readership of this article may get board to switch to the introduction section and again back to the result section. The active site domain of the protein is explained in Figure1, and a binding pocket of the complex protein-inhibitor is picturizing separately in Figures 2, 4, 6, and 8.

Comment#5-1:

I still think that a scheme of the inhibitors in the introduction section would increase the clarity of the study.

Answer:

Inhibitors' structure image was replaced in the introduction section (Fig2 and lines 62-69).

2) Answer #3:

For the "On the other hand, the inhibitor interacts such as independent fragments with different parts of the proteins. I think these findings are not very novel" we explained in the answer of Comment #16 of reviewe1. Also, we added a more detailed explanation of the Free energy landscape plots in the method section.

From the free energy landscapes, would be possible to predict a free energy of binding value of each system?

Answer:

For calculating binding free energy, there are some ways involving: alchemical transfer method (ATM), linear interaction energy (LIE), PDLD/S-LRA method, alchemical free energy perturbation (FEP), BAR method, and so on. But, We computed the free energy landscape (FEL) to capture inhibitors' lowest energy stable states during unbinding by gmx sham module (https://manual.gromacs.org/documentation/2018/onlinehelp/gmx-sham.html).

This module calculates free energy landscapes by computing the joint probability distribution from the two-dimensional plane constructed using two quantities (in our case, they were ligand RMSD (as the x-axis) and protein RMSD (as the y-axis)). Conformations sampled during the simulation were projected on this two-dimensional plane, and the number of points occupied by each cell was counted. The grid cell containing the maximum number of points is then assigned as the reference cell, with a free energy value of zero. Free energies for all the other cells were assigned with respect to this reference cell using the following equation:

ΔG(x,y) = -KbT ln P(x,y)/Pmin

P(x,y) estimates the probability density function obtained from a histogram of MD data, and Pmin is the maximum of the probability density function. Kb is the Boltzmann constant, and T is the temperature corresponding to each simulation.

3): The initial Molecular Dynamic Simulations should be longer to obtain an equilibrated configuration of the proposed models. Before SuMDs, the system should be very well equilibrated. You can see many studies about 3CL enzyme in which classical MD of microsecond timescales is performed before to study binding or reaction processes. The initial protocol of this work does not allow, perhaps, to explore important configurations of the enzyme that should be considered in the study. Authors should explore this way carefully for future studies.

Answer:

Thanks for your valuable advice, and we will consider them for our future research.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

17 Jan 2022

Reconstruction of the unbinding pathways of noncovalent SARS-CoV and SARS-CoV-2 3CLpro inhibitors using Unbiased Molecular Dynamics

PONE-D-21-32120R2

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,

Chandrabose Selvaraj, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

**********

4. 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

28 Jan 2022

PONE-D-21-32120R2

Reconstruction of the unbinding pathways of noncovalent SARS-CoV and SARS-CoV-2 3CLpro 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. Chandrabose Selvaraj

Academic Editor

PLOS ONE