The crystal structures of FP-2 and FP-3 were obtained from the protein databank (PDB) with the IDs 3BPF and 3BPM, respectively. Co-crystallized water molecules and other heteroatoms were removed, and Hydrogen atoms were added to the structures using the H++ server (Gordon et al., 2005). The 3D coordinates of the compounds QOD and ICD were collected from the MolPort database in SDF format, protonated, and converted to PDB format using the Schrodinger Maestro suite (Schrödinger, 2022). The cleaned-protonated proteins and ligands structures were then used to generate topology in the amber forcefield using AmberTools22 (Case et al., 2017). In each case, the ligand (i.e., QOD or ICD) is placed in a random position away from the protein (FP-2 or FP-3) in a cubic simulation box (Figure 2B) and solvated in the TIP3P water model (Price and Brooks, 2004). The systems were further neutralized with 150 mM NaCl, and ligand parameters were defined using GAFF forcefield (Wang et al., 2004).

All MD simulations were performed using GROMACS-2021.3 (Berendsen et al., 1995). Accordingly, each system was energy minimized for 5000 steps with the steepest descent algorithm followed by equilibration in the NVT ensemble for 100 ps to a temperature of 298 K using the velocity rescale temperature coupling (Bussi et al., 2007). The systems were then equilibrated in an NPT ensemble for 2 ns to a pressure of 1 atmosphere, using the Parrinello-Rahman barostat for pressure coupling (Parrinello and Rahman, 1998). This was followed by 1 𝛍s long final production mdrun in NPT ensemble, saving snapshots after every 10 ps. The complete simulation input parameters are described in the Supplementary Material.

To establish the stability of the systems, the trajectories generated from the 1 𝛍s simulations were used to calculate the root-mean-square deviation (RMSD) of the protein’s backbone atoms. Unless told otherwise, the last 500 ns of the trajectories were used for other downstream analyses, such as the root-mean-square fluctuations (RMSF). Clustering analysis was performed to classify the different binding orientations of the ligands. The clustering was done using the GROMACS gmx cluster module, and the GROMOS method was chosen as the classification method (Daura et al., 1998). Clusters are separated by a 0.25 nm difference in RMSD. The calculation was performed from a snapshot every 10 ps.

The binding energies of the interaction between QOD, ICD, and FP-2 or FP-3 were calculated from the last 500 ns of the MD simulation trajectories using the molecular mechanics generalized Born and surface area continuum solvation approach. In this method, the binding free energy is given by: Δ G b i n d = Δ H − T Δ S ⋍ Δ E M M + Δ G s o l − T Δ S (1) Δ E M M = Δ E int e r n a l + Δ E e l e c t r o s t a t i c + Δ E v d w (2) Δ G s o l = Δ G G B + Δ G S A (3)

(Hou et al., 2011)

The term ΔE_MM represents energy due to change in the gas phase, which is defined by internal (bond, angle, and dihedral energies, ΔE_internal), electrostatic (ΔE_{electrostatic}), and van der Waals energies (ΔE_vdw), while ΔG_sol and TΔS account for the free energy of solvation and entropy due to binding-induced conformational changes, respectively. The solvation free energy is a term that is derived from the sum of electrostatic (ΔG_GB, estimated from the Generalised Born model) and non-electrostatic (ΔG_SA) solvation energy components, which is calculated using the solvent accessible surface area (Genheden and Ryde, 2015). The binding free energy calculation was performed using gmx_MMPBSA software, and the input parameters are described in the Supplementary Material (Valdés-Tresanco et al., 2021).

The root-mean-square deviation (RMSD) is often used to measure global conformational changes in macromolecular structures and has become an increasingly popular method for assessing the convergence of molecular simulations. The RMSDs of the backbone atoms of our systems were calculated with reference to the energy-minimized structure, which is very close to the crystal structures. Although the systems were relatively stable even before the first 200 ns of the simulation, we observed some distortions around 500 ns simulation time, for example, in the FP-2-ICD complex (Figure 3A). This indicated that the system has not fully stabilized even at 500 ns, and therefore, those time frames should not be considered for analyses. In comparison, all systems with FP3 were remarkably quite stable throughout the simulation time, suggesting less dynamic interaction of FP3 with the ligands. Nevertheless, we deemed the first 500 ns as an extension of equilibration and only considered the last 500 ns for further structural analyses. The local protein flexibility can be followed using per-residue root-mean-square fluctuation (RMSF) and can especially be informative in describing ligand-induced flexibility. From Figure 3B, we can observe dramatic fluctuations around residues 107–120, 222–226 in the FP2-QOD system and residues 59–62, 79–83 in the FP-2-ICD complex, compared to the apo FP-2. These residues constitute part of the α4/β2, β6-β7, α2-α3, α3/α4 loops (respectively) in FP-2 (Hogg et al., 2006), and therefore fluctuations are expected. Other observed regions of flexibility include residue 178 to 183 in FP2-QOD, which constitute a helical part of the protein. Therefore, the ligand, through its interaction with the protein, could also contribute to the more pronounced fluctuations by influencing the global fluctuation of the system. Similarly, a higher fluctuation of the loops formed by residues 88–91, 105–130, and 232–242 (belonging to α3/α4 and β6-β7 loops) was observed in the FP-3-ICD complex, which could be due to the same aforementioned reason (Figure 3C). On the contrary, the RMSF profile of the FP3-QOD system correlates very well with apo-FP3 RMSF, suggesting that QOD may not have a significant effect on the local flexibility of the protein.

To understand the binding mechanism of the ligands to FP-2, RMSD-based clustering was applied to the generated trajectories to classify the different conformations visited during the simulation. The classification produced 278 and 188 clusters for FP2-QOD and FP-2-ICD, respectively. In each case, the largest clusters (31.13% and 41.96% of the total frames for FP2-QOD and FP2-ICD) happened to be the only cluster in which the ligand is bound to the target. Interestingly, both ligands occupy the same binding pocket at the active site of FP-2 (Figure 4A) and are coordinated by almost identical FP-2 residues (Figures 4B, C). QOD interacts with Q36, N38, A157, W206, Q209, and W210, with hydrogen bonds formed by Q209 and W206 with one of the oxygen atoms of the dioxanyl ring and the carbonyl oxygen of the ethanediamide ligand backbone, being the major stabilizing interactions (Figure 4B). Similarly, ICD interacts with Q36, A157, W206, and W210, in addition to D35 and K37, both of which are responsible for hydrogen bonding with the formamide carbonyl oxygen of the ligand. An additional hydrogen bond is formed between the side chain carboxyl oxygen of D35 and the nitrogen atom of the indoyl ring of the ligand (Figure 4C). FP-2 is a member of the papain-like C1 cysteine proteases family that is characterized by catalytic quads of cysteine, histidine and asparagine, and glutamine (Martynov et al., 2015), represented by Q36, C42, H174 and N204 in FP-2 (Ettari et al., 2010). Our predicted interaction models suggest that access to these catalytic residues by the natural substrate is prevented by the ligands.

Similarly, the trajectories from the simulation of the FP3-QOD and FP3-ICD complexes were clustered to classify the different conformations sampled. Just like in FP-2, the largest clusters (67.54% and 26.62% of the total frames for FP3-QOD and FP3-ICD, respectively) happened to be the only clusters where both ligands are bound to FP-3. Unlike in FP-2, the ligands are not fully embedded in the FP-3 active site and do not assume the same binding conformation but are docked very close to the critical residues, with the indonyl ring of ICD directly blocking the catalytic N213 from the top (Figure 5A). The interaction of FP-3 with QOD is maintained by a hydrogen bond between the side chain nitrogen atom of N86 and amide-oxygen of the ethanediamide backbone, as well as two hydrogen bonds formed between the backbone oxygen atom of Y90 and both N′-nitrogen and N-oxygen atom of the ligand ethanediamide backbone (Figure 5B). Other FP-3 residues that coordinate QOD include catalytic C51 and H183, Y93, I94, N96, S158, A184, and E243. The catalytic residues of the protein are not directly blocked in this interaction mode. However, the ligand’s methyl-tetrahydroquinolinyl ring is flexible and could move back and forth or even sideways to cause steric hindrance and prevent access to C51 and H183 for enzyme catalysis.

On the other hand, ICD forms only a single hydrogen bond with the side chain oxygen atom of D163 and the nitrogen atom of the ligand indonyl substituent. The amino acids A46, A166, N182, catalytic H183, W215, and W219 constituted the remaining residues interacting with ICD (Figure 5C). Notably, the ligand’s biphenyl rings appeared not to contribute significantly to this interaction, given that the rings protruded outside and were almost completely excluded from the protein surface (Figures 5A, C).

The calculated binding energies of QOD and ICD against FP-2 returned comparable values (Table 1). This is not surprising, given that both ligands occupy the same binding pocket at the receptor’s active site and are coordinated by almost identical residues (Figure 4). This is also consistent with the experimental IC₅₀ values, which are 14.71 μM and 12.16 μM against FP-2, for QOD and ICD, respectively (Rana et al., 2020). The difference in the predicted binding energies of the ligands against FP-3 is surprisingly large (Table 1), considering their comparable IC50 values of 9.98 μM and 8.56 μM for QOD and ICD, respectively. However, this observation agrees, to some extent, with the predicted interactions since QOD exhibited more robust interaction with FP-3 in terms of the number of hydrogen bonds and the coordinating residues (Figure 5).