Structure retrieval and Modeling of Peptide (ATLQAIAS)
The crystal structure of the 3CLpro complex was retrieved from the protein databank ( http://www.rcsb.org/) and assessed for Structural deformities. The peptide sequence (ATLQAIAS) was obtained from the previous literature and modeled using PEPFOLD3(https://bioserv.rpbs.univ-paris-diderot.fr/services/PEP-FOLD3/), which is an online web server for peptides modeling. Previous computational and experimental studies reported that this peptide, along with others, remained active against the SARS virus and potentially inhibited the virus replication[20, 21].
To test the activity of the (ATLQAIAS) peptide against the SARs-CoV-2, molecular docking of the peptide and 3CLpro from SARs-CoV-2 was carried out to understand the interaction pattern. For this purpose, ZDOCK (http://zdock.umassmed.edu/) was used with a blind docking option to increase and test the reliability of the server. The top ten docked conformations were retrieved and analyzed for the interaction with the key residues such as His41 and Cys145, which form a catalytic dyad and do the proteolysis.
Interface analysis and Peptide library construction
Understanding the binding interface of two interacting protein or peptide is vital for designing small peptide inhibitors. The use of small peptides derived from the native peptide has been widely used in the designing of short peptide derivatives. A similar strategy to manipulate the interface residues of the reference-native-peptide-3CLpro and to design decoy peptides against SARs-CoV-2 using the machine learning protocol implemented in MOE. The initial docking of the peptide (ATLQAIAS) with the 3CLpro confirmed that it also blocks the catalytic dyad and form supplementary interaction with the key active site residues. After this confirmation, information regarding the key residues in the interface is important for mutants modeling. For this purpose, the Alanine scanning strategy was exercised to calculate the impact of each residue in the interaction with the 3CLpro. To grasp this strategy, the dAffinity and dStability parameters were monitored in the ASM (alanine scanning module) of MOE. The dAffinity and dStability values show the relative change in binding energy when a particular amino acid is changed into alanine. These scores reveal important information regarding the importance of a particular residue, after the confirmation of key residues selection of the hotspot residues for residue scan to be replaced by the 19 amino acids. For this purpose, the residue scan module of MOE was used using the Unary Quadratic Optimization (UQO) parameter under the LowModeMD ensemble. This residue scan generated a database of mutant peptides with their respective scores. The detailed mechanism of this alanine scanning mutagenesis and residue scan approaches has been discussed previously.
Molecular Dynamics Simulation
A comprehensive protein dynamics technique was used to evaluate the dynamic features of the 3CLpro and the designed peptides via the residue scan method. Rationally developed peptides with fierce dAffinity and dStability scores were obtained using the Amber18 package and subjected to molecular dynamics simulation (MDS) with the ff14SB force field . Using the TIP3P water model with box dimension 10.0Å, each structure was correctly solved and neutralized by incorporating Na+ counter ions. Two steps of energy minimization were performed. For the first energy minimization, 6000steps, while for the second 3000 steps of conjugate gradient minimization, was completed. Following the minimization, each system was heated using default parameters (300 K for 200 ps). Weak restraint was used for density equilibration for 2ns, while the whole system at a constant pressure for 2ns. A 100ns MD under constant pressure was performed. For the temperature control, Langevin thermostat (1 atm, 300 K) was used . Particle Mesh Ewald (PME) algorithm was used to compute long-range interactions [27, 28]. The cutoff distances were set to 10Å. For the covalent bonds involving hydrogen, the SHAKE algorithm was used. GPU accelerated simulation using (PMEMD.CUDA) was used for all the processes.
Post simulation analyses such as RMSD, RMSF, and radius of gyration (Rg) were calculated by using CPPTRAJ and PTRAJ modules of Amber .
Binding affinity calculations
To estimate the binding free energy of the designed peptides towards the 3CLpro MMPBSA.PY script was used (Miller et al., 2012). Using 5000 frames from the 100ns trajectory, the free energy of binding was estimated by using the following equation:
∆𝐺𝑏𝑖𝑛𝑑 = ∆𝐺𝑐𝑜𝑚𝑝𝑙𝑒𝑥 − [∆𝐺𝑟𝑒𝑐𝑒𝑝𝑡𝑜𝑟 + ∆𝐺𝑙𝑖𝑔𝑎𝑛𝑑]
In the above equation, the ∆Gbind represent the total BFE while the other components in the equation represent BFE of the complex, the protein, and ligand molecules. Each term in the binding free energy was estimated using the following equation:
𝐺 = 𝐺𝑏𝑜𝑛𝑑 +𝐺𝑒𝑙𝑒 +𝐺𝑣𝑑𝑊 +𝐺𝑝𝑜𝑙 +𝐺𝑛𝑝𝑜𝑙 −𝑇𝑆
The above equation represents the polar, nonpolar, electrostatic, solvent-accessible surface area SASA, and van der Waals interactions, respectively.
Per-residues Energy Decomposition analysis
Per-residue energy decomposition is the best approach to understand the energetic contribution at the residues level significantly. Herein, to understand the impact of each substitution on the binding of the designed peptides, we used per-residue energy decomposition analysis.
Clustering of MD trajectories using PCA
To comprehend the motion of MD trajectories, an unsupervised learning method known as Principal Component Analysis (PCA)[31, 32]was performed to acquire knowledge regarding the internal motion of the system. For this purpose, an Amber module known as CPPTRAJ was used. The spatial covariance matrix was determined for eigenvector and their atomic co-ordinates. Using the orthogonal coordinate transformation, a diagonal matrix of eigenvalues was generated. Based on the eigenvectors and eigenvalues, the principal components were extracted. Using these PCs, the dominant motions during the simulation were plotted[33, 34].