## 2 Materials and Methods  
In this work, we have used a combination of flexible molecular docking and Density Functional Theory (DFT) calculations. Docking methodologies were used to reasonably estimate the adduct between esbelen and \(\mathrm{M}^{\mathrm{pro}}\). The mechanism was investigated employing the so-called cluster approach [21], which is the fully quantum investigation of the enzyme reactivity using only selected residues that reproduce the chemical features of the catalytic pocket. Such a methodology was largely applied to the reactivity of highly different enzymes, spanning oxidoreductase [22], lyase [23] and serine [24] or metalloprotease [25] enzymes. While the inclusion of the enzymatic environment at a Molecular Mechanics level of theory (i.e., employing a QM/MM methodology such as the ONiOM scheme) [26]) is expected to have an impact on the reactions' energetics, [PERSON] and coworkers [27] showed that the surrounding effect is rather weak when the catalytic pocket is exposed to the solvent, as it is for SARS-CoV-2 \(\mathrm{M}^{\mathrm{pro}}\). Thus, the environment is expected to affect the energetics, but is unlikely to alter the mechanistic features of the reaction. Molecular docking was carried out to simulate the binding pose of EbSe with \(\mathrm{M}^{\mathrm{pro}}\) from SARS-CoV-2 (PDB ID 6 LU7 [6]). AutoDock Vina software [28], which has a high accuracy for binding mode predictions, was used [29]. For the \(\mathrm{M}^{\mathrm{pro}}\) structure preparation, the waters, ions, ligands, and other molecules were removed, while the hydrogen atoms were added using the CHIMERA program, followed by 100 steps of energy minimization (amberff99 SB) [30]. The catalytic dyad (Cys145 and His41) was considered neutral, as highlighted in previous studies [31; 32]. An adduct between EbSe and cysteine (EbSe-Cys) was used in the docking simulations to mimic the putative metabolites of EbSe [33; 34; 35]. The tridimensional model of the ligands (EbSe and EbSe-Cys) were created with Avogadro and MOPAC (PM6 method) [36; 37]. The files were prepared for the docking, using the AutoDock Tools 4.2 [38], with the flexible ligands. To consider the protein induced-fit effect and to improve the interactions between ligands and \(\mathrm{M}^{\mathrm{pro}}\), the flexible-flexible docking method was applied, where EbSe, EbSe-Cys, and the side chain of His41, Met49, Asn142, Cys145, Met165, Glu166, and Gln189 residues (from the active site) were considered flexible during the simulations. The grid box was centered on the coordinates x = \(-\)14.04, y = 17.44, and z = 66.22 (size = 25 \(\times\) 25 \(\times\) 25 A), and an exhaustiveness of 50 was used. The complex ligand-receptor with the most favorable binding free energy, and the best Se\(\cdots\)S(Cys145) orientation, was selected as a model of the binding pose and used to build the \(\mathrm{M}^{\mathrm{pro}}\) pocket cluster for the DFT study. After the docking, taking into account the closest \(\mathrm{M}^{\mathrm{pro}}\) residues (at 4.0 A) from EbSe and EbSeS, His41, Met49, Cys145, and Met165 were selected and removed from the \(\mathrm{M}^{\mathrm{pro}}\) pocket, and the CH\({}_{3}\)CO and CH\({}_{3}\)NH groups were added to the N- and C-terminal regions, respectively, to mimic the backbone peptide bonds. In fact, the His41, Met49, Cys145, and Met165 residues were involved in many interactions with small molecules [39]. As a general model of esbelen thiolates, the carboxyl and amino moieties of EbSe-Cys were replaced with H atoms to create the ethanthiolate (EbSeS). All density functional theory calculations were performed, employing Gaussian16 [40]. The level of theory employed was based on previous mechanistic investigations in an enzyme whose catalysis is based on reactive Cys or Sec residues [41; 42]. Thus, the B3 LYP hybrid functional [43; 44], in combination with Grimme D3 dispersion correction and the Becke-Johnson damping function [45; 46] was used.