2.  Materials and methods
In general, molecular modelling is implemented as an essential tool for the prediction of drug–macromolecule interaction. This technique helps to enhance the success rate of an experiment and cuts down the experimental cost. Hence, the molecular docking study can help to analyse the possible binding pose of a small molecule on the active site of a macromolecule. Here we used molecular docking to screen some biologically active spice molecules with the SARS-CoV-2 RBD Spro and SARS-Cov-2 Mpro. The molecule with the highest binding affinity to RBD Spro and Mpro was subjected to MD simulation for further validation.

2.1.  Drug-likeness properties of the small molecules
The property of the small molecules for drug-likeness was estimated using the Lipinski’s rule (Lipinski et al., 2001). This rule works on five parameters viz. no more than five hydrogen bond donors, no more than 10 hydrogen bond acceptors, molecular mass <500 Da, and the octanol-water partition coefficient, i.e. log P should not exceed 5. The Lipinski’s parameters were obtained by using the SwissADME server (www.swissadme.ch/index.php) (Daina et al., 2017).

2.2.  Structure preparation of the proteins and ligands
The crystal structures of the SARS-CoV-2 Spro (PDB ID: 6M0J) and SARS-CoV-2 Mpro (PDB ID: 6Y84) were obtained from the RCSB protein data bank. All the non-standard residues, including water, were removed from the PDB file using Chimera (Pettersen et al., 2004). The RBD domain was obtained from the PDB file of SARS-Cov-2 Spro. The RBD Spro residue sequence number before and after removing ACE2 is presented in Supplementary Figure S1. The 3D conformers of the ligands were obtained from the PubChem and optimized using the steepest descent and conjugate gradient steps with General Amber Force Field (GAFF) (Wang et al., 2004) in Chimera (Pettersen et al., 2004).
Figure 1.  Predicted lipophilicity (Log P) values of the spice molecules obtained from different calculation models.

2.3.  Molecular docking study
The prepared structures of the protein and ligand were subjected to molecular docking analysis using AutoDock Vina (Trott & Olson, 2010). AutoDock Vina is the newest member of the AutoDock family that has improved speed and accuracy. It uses a hybrid scoring function and a quasi-Newtonian optimization algorithm to find the lowest energy confirmations within the search space. A grid box of 40 Å × 65 Å × 70 Å was built with the centre of the box at (11.98, 0.60, 4.79) for the SARS-CoV-2 Mpro. A grid box of size 30 Å × 45 Å × 30 Å with centre at (−36.51, 30.69, 5.48) was prepared for the SARS-CoV-2 RBD Spro. The exhaustiveness of search was set at 20 and 8 for the SARS-CoV-2 Mpro and the SARS-CoV-2 RBD Spro, respectively, to compensate for the larger box volume and reliable results. The docked poses were ranked as per their binding affinities at the end of the docking run. The ligand interactions of the best-docked poses at the active sites of the macromolecule were extracted using PyMol (Schrödinger LLC, 2017). The ligand interactions were analysed using the 2D interaction plot in the Discovery Studio Visualizer (2005). The Coulombic electrostatic potential surface was determined with the help of the APBS plugin available in PyMol (Schrödinger LLC, 2017).

2.4.  Molecular dynamics (MD) simulation study
To verify the stability of the complex and interaction dynamics of SARS-CoV-2 RBD Spro and SARS-CoV-2 Mpro, we performed the MD simulation study of the two complexes with the highest docking score (Spro-Piperine and Mpro-Piperine), using GROMACS-5.1.5 (Abraham et al., 2015). The CHARMM36 forcefield (Huang & MacKerell, 2013) was used for the simulation of the systems. The topology parameters for the ligand molecule were obtained from CGenFF (Vanommeslaeghe & MacKerell, 2012). A dodecahedron simulation box filled with TIP3P water model (Price & Brooks, 2004) was prepared. Counter ions were added to maintain the electrical neutrality of the system. The systems were kept at a buffer concentration of 0.15 M. Then, the build systems were energy minimized with 50,000 steps using the steepest descent algorithm. Then the systems were equilibrated under NVT and NPT ensembles for 100 ps at 300 K temperature and 1 atm pressure before the production run. After the equilibration, a production run of 100 ns was incorporated under the NPT ensemble. For long-range electrostatic interaction, particle mesh Ewald (PME) (Darden et al., 1993) and for van der Waals interactions, the force-switching scheme was used. Besides, for temperature and pressure coupling, the Berendsen thermostat (Berendsen et al., 1984) with velocity rescaling and Parrinello-Rahman barostat (Parrinello & Rahman, 1981) with isotropic rescaling were used, respectively. The simulation time step was set to 2 fs, and the trajectories were recorded at every 10 ps. The simulation data were analysed by analysing the root mean square deviation (RMSD), root mean square fluctuation (RMSF), number of hydrogen bonds and radius of gyration (Rg) using Gromacs analysis tools. The principal component analysis (PCA) was performed using the g_covar, g_anaeig and g_sham tools of Gromacs. The data were exported to origin 9.0 and plotted for further analysis purposes.

2.5.  Mm/PBSA binding free energy calculation
The method of calculation of binding free energy from MD trajectory snapshots using the molecular mechanics Poisson–Boltzmann surface area method is widely used. The binding free energy of the systems was estimated by extracting the snaps from the last 20 ns of the MD simulation using g_mmpbsa tool of Gromacs (Baker et al., 2001; Kumari et al., 2014). The binding free energy takes the contribution from vacuum potential energy, polar solvation energy and non-polar solvation energy. The binding free energy can be represented as (1) ΔGbind=Gcomplex−(Gprotein+Gligand) where Gcomplex, Gprotein and Gligand are the total free energies of the complex, isolated protein and isolated ligand, respectively. The free energy of the individual terms was estimated by (2) Gx=EMM−TS+Gsolvation where x is the complex, protein or ligand, and TS represents the entropic contribution to free energy in a vacuum with T and S as temperature and entropy. The average molecular mechanics potential and solvation free energies were calculated by using Equations (3) and (4) (3) EMM=Ebonded+Enonbonded= Ebonded−(Eelec+Evdw) (4) Gsolvation=Gpolar+Gnonpolar where Ebonded takes the contribution from a bond, angle and dihedral terms and Enonbonded consists of electrostatic and van der Waals energy contributions. The solvation energy includes the polar and non-polar solvation energies from the Poisson–Boltzmann equation and solvent accessible surface area (SASA), respectively.