System Setup: Model structures of the protein-ligand complexes were downloaded from the Protein Data Bank. The PDB accession codes for the complexes studied include 6LZE, 12 7D1M, 15 7D3I, and 6XHM, 17 for 11a, GC373, MI-30, and PF-231, respectively. The covalent bond between the warhead of the inhibitor and active-site cysteine (Cys145) of the protein were broken prior to running the molecular dynamics (MD) simulations. The M pro amino acid residues were assigned protonation states corresponding to their default values at neutral pH. We have previously evaluated methods for the calculation of the pKa of titratable cysteine residues in proteins and have performed MD simulations to estimate the protonation state of the catalytic Cys145/His41 dyad of the M pro . Based on our results and the work of others, we treated the catalytic Cys145/His41 dyad of the M pro as neutral residues. We also modeled His41 residue in two tautomeric states (i.e., N-𝛿 and N-𝜀 protonated states) to explore the preferred protonation state. The protonation state of histidines in SARS-CoV-2 M pro , in particular His41, has been shown to impact ligand binding and protein structure dynamics. The GROMACS MD software 41 (version 2020.4) was used to perform the calculations. Model structures of the protein-ligand complexes were prepared in explicit solvent using the TIP3P water model. The AMBER99SB-ILDNP and GAFF 44 force fields were used to generate parameters for the protein and ligands, respectively. Ligand charges were parameterized using AM1-BCC. The model systems were solvated in a dodecahedral box with a maximum cutoff distance of 12 Å from the edge of the box. The simulation cell was charge neutralized with an appropriate number of Na + ions. Following this, energy minimization and 2-ns equilibration runs (using both NVT and NpT ensembles) were performed on the system. This was subsequently followed by 10-ns simulation of the complex. The MD simulations were performed to determine the preferred protonation state of His41 in the model complexes prior to running the binding free energy simulations.