3.2.1. Molecular Docking

The docking study was carried out against both the receptor-binding domain (RBD) of SARS CoV-2s spike protein (S-protein) and human angiotensin-converting enzyme-2 (ACE2) (PDB codes: 6lZG and 1R42, respectively) [113,114]. An Autodock Vina docking machine [105] was used for docking experiments. To determine the best binding site for docking, we docked heparin (i.e., sulphated poly saccharide analogue) against both S-protein and ACE2 using the ClusPro server [14] by uploading each crystal structure to the website and choosing heparin as a ligand. The retrieved docking poses were then arranged from the highest to the lowest scores. We selected the docking pose of the highest score and subjected it to a short molecular dynamic simulation (25 ns) experiment to study the stability of heparin in this pose. Heparin showed acceptable stability in this orientation, where it deviated from the initial pose by an average RMSD of 2.1 Å (Figure 6); hence, the binding site of heparin in this pose was used for docking in the subsequent docking experiments. Heparin docking was carried out using its tetrasaccharide form (i.e., the default form in the ClusPro software). Accordingly, we also used the tetrasaccharide form of each MSP for docking experiments to reduce the computational cost required for docking and molecular dynamic simulations. All results were visualized using PyMol software [115].

#### 3.2.2. Molecular Dynamic Simulation

MD simulation experiments were carried out using the MDS machine in Maestro software, Desmond v. 2. [116–118], using its default force field (i.e., OPLS-AA). Protein– ligand systems were constructed using the System Builder function. Thereafter, these systems were embedded in an orthorhombic box consisting of TIP3P water and 0.15 M Na+ and Cl- ions (the default dimensions were used). Subsequently, the prepared systems were energy-minimized and equilibrated for 10 ns. Ligand parameterization was carried out during the system building step according to the OPLS force field. For MD simulations carried out using NAMD software, the parameters and topologies of the ligands were computed using the Charmm36 force field. The online software Ligand Reader and Modeler (http://www.charmm-gui.org/?doc=input/ligandrm, accessed on 15 May 2021) [119] and the VMD plugin Force Field Toolkit (ffTK) [120] were utilized for this regard. The generated parameters and topology files were then loaded into VMD to read the protein– ligand complexes without errors, then the simulation step was performed. MD simulations were run for 50 ns at 310 K in the NPT ensemble with the Nose–Hoover thermostat and Martyna-Tobias-Klein barostat using anisotropic coupling. We selected the best binding poses for each compound as starting co-ordinates to investigate their binding stability and mode of interaction.

#### 3.2.3. Binding Free Energy Calculations

Binding free energy calculations ( Δ *G*) were performed using the free energy perturbation (FEP) method [121]. Briefly, this method estimates the binding free energy (i.e., Δ Gbinding) according to the following equation: Δ *<sup>G</sup>*binding = Δ *<sup>G</sup>*Complex − Δ *<sup>G</sup>*Ligand. These estimations were derived from separate simulations (NAMD software was used for these experiments). All input files required for simulation by NAMD were papered using the online website Charmm-GUI (https://charmm-gui.org/?doc=input/afes.abinding, accessed on 30 April 2021). Subsequently, we loaded these files into NAMD in order to produce the required simulations. The FEP function in NAMD was used to accomplish this experiment. The equilibration step was achieved in the NPT ensemble at 300 K and 1 atm (1.01325 bar) with Langevin piston pressure (for "complex" and "ligand") in the presence of the TIP3P water model. Then, 10 ns FEP simulations were carried out for each ligand and the last

5 ns of the free energy values was measured for the final free energy estimation [121]. All resulting trajectories were visualized and analyzed by VMD software.
