Computational Design of Miniprotein Inhibitors Targeting SARS-CoV-2 Spike Protein.

The ongoing pandemic of COVID-19 caused by SARS-CoV-2 has become a global health problem. There is an urgent need to develop therapeutic drugs, effective therapies, and vaccines to prevent the spread of the virus. The virus first enters the host cell through the interaction between the receptor binding domain (RBD) of spike protein and the peptidase domain (PD) of the angiotensin-converting enzyme 2 (ACE2). Therefore, blocking the binding of RBD and ACE2 is a promising strategy to inhibit the invasion and infection of the virus in the host cell. In the study, we designed several miniprotein inhibitors against SARS-CoV-2 by single/double/triple-point mutant, based on the initial inhibitor LCB3. Molecular dynamics (MD) simulations and trajectory analysis were performed for an in-depth analysis of the structural stability, essential protein motions, and per-residue energy decomposition involved in the interaction of inhibitors with the RBD. The results showed that the inhibitors have adapted the protein RBD in the binding interface, thereby forming stable complexes. These inhibitors display low binding free energy in the MM/PBSA calculations, substantiating their strong interaction with RBD. Moreover, the binding affinity of the best miniprotein inhibitor, H6Y-M7L-L17F mutant, to RBD was ∼45 980 times (ΔG = RT ln Ki) higher than that of the initial inhibitor LCB3. Following H6Y-M7L-L17F mutant, the inhibitors with strong binding activity are successively H6Y-L17F, L17F, H6Y, and F30Y mutants. Our research proves that the miniprotein inhibitors can maintain their secondary structure and have a highly stable blocking (binding) effect on SARS-CoV-2. This study proposes novel miniprotein mutant inhibitors with enhanced binding to spike protein and provides potential guidance for the rational design of new SARS-CoV-2 spike protein inhibitors.

Introduction

The novel strain of coronavirus (CoVs), Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) was originally named the 2019-nCoV and has been confirmed as the pathogen of an ongoing global outbreak of respiratory disease that has caused serious damage to the world economy and global public health.[1−3] At present, seven human coronaviruses have been identified worldwide, including human coronavirus 229E (HCoV-229E), HKU1 (HCoV-HKU1), OC43 (HCoV-OC43), NL63 (HCoV-NL63), Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV), Middle East Respiratory Syndrome Associated Coronavirus (MERS-CoV), and SARS-CoV-2.[4]

Coronavirus are pathogens of the Coronaviridae family and have an impact on the health of humans and animals, especially the respiratory and gastrointestinal systems, with symptoms ranging from mild to fatal. Coronavirus are classified into the family Coronaviridae and the order Nidovirales.[5] Coronavirus are large, enveloped, positive single-stranded RNA virus, which can be further divided into α-, β-, γ-, and δ-CoVs. Of these four subtypes, α-CoVs and β-CoVs were identified to infect humans.[4] The COVID-19 pandemic was caused by SARS-CoV-2, which is an enveloped virus belonging to the Coronaviridae family, the Coronaviridae subfamily, and the β-Coronavirus genus.[6] Structurally, SARS-CoV-2 has the largest RNA genome among any other known viruses, ranging in size from 26 to 32 kb, with a particle diameter of 125 nm.[7] Its characteristics are nonsegmented, positive, and single-stranded RNA.[8] Meanwhile, its genome has 96.2% homology with bat CoV RaTG13, 79.6% homology with SARS-CoV, and ∼52% homology with MERS-CoV.[6,9] Despite a high degree of structural homology, SARS-CoV-2 is much more infectious and lethal than SARS-CoV and MERS-CoV. SARS-CoV-2 mainly causes severe acute and other complications in the kidney, spleen, heart, and brain, eventually leading to acute respiratory distress syndrome, severe hypoxemia, viral pneumonia, and gastrointestinal and neurological symptoms in patients with COVID-19.[10]

The genome of SARS-CoV-2 encodes four conserved structural proteins, including spike (S), membrane (M), envelope (E), and nucleocapsid (N).[11−13] The S protein on the surface of the virus particle is responsible for targeting and attaching the host cell and triggering the fusion of the virus with the host cell membrane, which is a key step in initiating infection and transferring viral RNA to the host cell. The M protein and E protein are associated with virus shape and assembly/germination, respectively. The N protein enters host cells, along with the genetic material of SARS-CoV-2, contributing to RNA transcription, replication, virus assembly, and release.[14,15] The genome of SARS-CoV-2 also encodes 16 nonstructural proteins (NSP1-NSP16). A typical coronavirus genome contains at least 6 open reading frames (ORFs), which are responsible for encoding these nonstructural proteins (NSPs).[16] ORF1a and ORF1b are translated into two large nonstructural viral polyproteins: pp1a (4382 amino acids) and pp1ab (7073 amino acids). They are processed by 3C-like protease (3CLpro) and papain-like protease (PLpro), respectively, and produce a variety of nonstructural proteins, including RNA-dependent RNA polymerase (RdRp) and helicases, to catalyze viral genome replication and protein synthesis.[17,19] In addition, the genome of SARS-CoV-2 also encodes nine putative accessory proteins. Most of the proteins encoded by SARS-CoV-2 are similar to SARS-CoV.[18,20] The S protein contains 1273 amino acids and consists of two subunits S1 (14–685 a.a.) and S2 (686–1273 a.a.), which are preceded by a short peptide (1–13 a.a.). The receptor binding domain (RBD) in the S1 subunit is responsible for direct binding to the peptidase domain (PD) of the angiotensin-converting enzyme 2 (ACE2) and mediating the entry of the virus into the host cell, while the S2 subunit is responsible for fusion with the host cell membrane.[21−31] The S1 subunit is composed of three domains: the N-terminal domain (NTD, 14–305 a.a.); the receptor binding domain (RBD, 319–541 a.a.); and the C-terminal domain (CTD), which contains two subdomains (SD1 and SD2).[10,32,33] The S2 subunit contains a fusion peptide (FP; 788–806 a.a.), which consists of the heptad repeat 1 (HR1; 912–984 a.a.), the heptad repeat 2 (HR2; 1163–1213 a.a.) and hydrophobic residues, a transmembrane domain (TM; 1213–1237 a.a.) and a cytoplasmic domain (CP; 1237–1273 a.a.).[10,33,34] In its natural form, the S protein exists as a trimer on the surface of the virus particle. The S1 and S2 subunits form the extracellular stalk and the bulbous crown.[10,30,33] Once the S protein binds to the receptor, the RBD of the S1 subunit will undergo a hingelike conformational movement, and the protein changes from the conformation before the fusion to the conformation after the fusion. RBD is a globular domain located on the distal surface of the spike protein. Two conformations are observed in the stable trimer. Specifically, in one conformation, an RBD is accessible to ACE2, while the other two are not, that is, the open conformation. In other conformation, all three RBDs are down, the ACE2 receptor is inaccessible, that is, the close conformation.[27,35] After binding to the ACE2 receptor, the spike protein trimer undergoes conformational changes.[36,37]

Over the past year or so, several classes of (bio)molecules, including nanobodies,[38−40] neutralizing antibodies,[41] soluble ACE2 protein derivatives,[42] peptides,[43] miniprotein,[44,45] and small molecules,[46] targeting the S protein have been designed and identified to prevent and block the novel coronavirus infection.[47,48] Among them, the current strategies for designing high-affinity miniprotein inhibitors are mainly as follows: first, most of the interaction regions with RBD from the protein binding domain of ACE2 were extracted to design miniprotein inhibitors or two successive α-helix bundles were extracted to design a miniprotein inhibitor by connecting them with an amino acid or carbon–carbon bond.[43,49,50] Second, de novo design of miniprotein inhibitors is independent of previous RBD-binding interactions.[44] Meanwhile, several studies have reported several bioactive molecules as well as drug repurposing molecules that will take multiple important proteins of SARS-CoV-2, including the main protease (Mpro), some nonstructural proteins, spike protein, and RNA-Dependent RNA Polymerase, as targets for which inhibitor screening was performed based on the in-silico techniques.[51−58] However, the studies using the point mutation and molecular dynamics (MD) simulations to screen high-affinity miniprotein inhibitors are limited.

In our study, several miniprotein inhibitors against SARS-CoV-2 by the single/double/triple-point mutant were designed on the basis of the recently announced crystal structure (PDB ID 7JZM).[44] The MD simulations and binding free-energy calculations, along with a variety of analysis methods, were performed to clarify the underlying mechanism of the interaction between RBD and inhibitory miniproteins, and investigate the effects of mutations on inhibitory activity of the designed inhibitors. Therefore, some new RBD inhibitors are put forward and the study provides potential guidance for the rational design of SARS-CoV-2 spike protein inhibitors.

Materials and Methods

System Preparation

In this study, the LCB3-RBD (PDB ID 7JZM)[44] complex structure was retrieved from the Protein Data Bank (www.rcsb.org), which was used as the template for a series of mutations. Initially, the mutants, H6Y, F30Y, L17F, H6Y-L17F, H6Y-M7L-L17F, were generated using Accelry’s Discovery Studio (DS).[59] The protonation states of the ionizable residues were predicted using the H++ web server.[60]

On the basis of the initial inhibitor LCB3, 64 residues were first mutated to alanine by using DS software. The mutation energy and the effect of mutation obtained from Calculate Mutation Energy (Binding) protocol in the DS are shown in Table S1 in the Supporting Information. Since the mutation of the 8 amino acids in the table leads to a decrease in affinity, these amino acids may be the key amino acids for the interaction of the inhibitor with the receptor RBD protein. As we all know, there are 20 common amino acids in nature. Hence, each amino acid can be turned into its 19 mutants. To improve the inhibitory activity of the inhibitor LCB3, these 8 amino acids were subjected to single-point mutations into an additional 19 amino acids. Theoretically, 152 single point mutants have been obtained and the mutation energies are shown in Table S2 in the Supporting Information. The mutants with His6 mutated to Tyr6 and Leu17 mutated to Phe17 have relatively low mutation energies of −1.39 and −1.08 kcal/mol, respectively. Furthermore, the mutant F30Y was selected because of the relatively low energy value of the binding free energy of the MD simulation, relative to the initial system. To further improve the inhibitory activity, double and triple point mutations were performed, based on the results of single-point mutations. According to Table S3 in the Supporting Information, the double-point mutant H6Y-L17F (−2.31 kcal/mol) and the triple point mutant H6Y-M7L-L17F (−3.11 kcal/mol) were also selected for further studies to determine the best inhibitor.

Molecular Dynamics (MD) Simulations

All MD simulations were performed using the GROMACS 4.5.3 package.[61] The Amber99SB-ILDN force field[62] was used to provide the parameters of the bonded and nonbonded interaction. The topologies and coordinates files of protein complexes were generated by the pdb 2gmx module of the GROMACS software. Under periodic boundary conditions, the three-site water model TIP3P[63] was used to solvate the protein complexes, forming an 87.45 Å × 87.45 Å × 87.45 Å cubic box with the nearest distance of 10 Å between protein complexes and the edge of the box. The counterions were added to neutralize the system, and the salt concentration of 0.15 M was set to simulate a physiological environment.

The energy minimization of all systems was employed by using the steepest descent algorithm. Then, the minimized system was heated gradually to 300 K for 100 ps of simulation by the v-rescale temperature coupling method[64] with a coupling constant of 0.1 ps in the canonical ensemble (NVT). Afterward, 1 ns simulation was performed to maintain the system pressure of 1 atm by the Parrinello–Rahman barostat[65] with a coupling constant of 2 ps and compressibility of 4.5 × 10–5 bar–1 in the constant pressure ensemble (NPT). Finally, each system was performed to the production simulation of 2 × 400 ns without any position constraints, in order to verify the reliability of the simulation results. During the simulations, the particle mesh Ewald (PME)[66,67] was employed to estimate the contributions of long-range electrostatic interaction. Furthermore, van der Waals interactions were computed with a cutoff distance of 10 Å. The LINCS algorithm[68] was applied to constrain all bonds concluding hydrogen atoms. A time step of 2 fs was set and the trajectory coordinates were saved every 10 ps. The Visual Molecular Dynamics (VMD) program[69] was applied to visualize the MD trajectories.

Interaction Analysis

In VMD, the thresholds of the D–A distance and D–H–A angle for hydrogen bonds (H-bonds) were set to 3.5 Å and 150°, respectively. The cutoff distance between the N atoms of the basic residue and the O atoms of the acidic residue was set to 3.5 Å for the salt bridge analysis.

Principal Component Analysis (PCA)

Principal component analysis (PCA)[70,71] is an effective method to reduce the high-dimensional data into the major components to reflect the main internal movement, and then recognize the dominant conformational change of a molecule during MD simulation. Each frame of the MD simulation trajectory was first fitted to an average structure to remove the translation and rotational motion of all Cα atoms. Then, a 3N × 3N covariance matrix was established with Cartesian coordinates, and the eigenvectors and eigenvalues were obtained by diagonalizing the matrix. The eigenvalues represent the amplitude of the motion and the eigenvectors represent the direction. In this work, porcupine plots are used to display the dominant motion of the inhibitor based on the first principal component.

Dynamic Cross-Correlation Matrix (DCCM)

Dynamic cross-correlation matrix (DCCM) was implemented to evaluate the dynamic correlation differences of several systems. In the biomolecular system, the correlated motion of Cα atoms is crucial in identifying binding regions. DCCM gives the pairwise correlations of the movement of amino acid residues in proteins during MD simulations and was calculated via the Bio3d package in R studio,[72] based on the simulation trajectories. In this work, the cross-correlation matrix between residues i and j was computed using the stable conformations of MD trajectories. The dynamics correlation of the RBD-inhibitor complexes was determined by analyzing the correlated movement of residues. The correlative motion degree of two residues (two atoms) can be expressed by the coefficient of cross-correlation C, which is defined aswhere i(j) represents the ith (jth) atom (or residue), Δr (Δr) corresponds to the displacement of the ith (jth) atom (or residue) mean position, and angle brackets indicate an ensemble average. The C value is in the range from −1 to 1. The positive value describes the positively correlated motion (the same direction), while the negative value means the anticorrelated motion (the opposite direction) within the period of simulations.

MM/PBSA Calculations

The calculation of the binding free energy is a vital thermodynamics method that can provide important details of the receptor–ligand interaction.[73] The molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) method was used to evaluate the binding energy of the initial inhibitor LCB3 and RBD and the effect of the mutant inhibitors on the binding energy of RBD from an energy perspective. The MM/PBSA assessment was performed using the g_mmpbsa tool[74] in the scripts-based APBS program. Based on the MM/PBSA method, the prediction of binding free energy between the RBD and miniprotein inhibitors can be achieved as described in the following equation:where ΔGbind represents the binding free energy between the RBD of spike protein and miniprotein inhibitors. The terms Gcomplex, GRBD, and Ginhibitor represent the free energy of the RBD–inhibitor complexes, RBD, and inhibitors, respectively. Here, the term ΔGbind is calculated using the following equations:where ΔGbind is composed of the energy of the gas phase (ΔEgas), solvation free energy (ΔGsol), and the contribution of entropy (−TΔS). The entropy is not considered due to the large computational cost and low prediction accuracy. ΔEgas is the result of the sum of the van der Waals interaction energy (ΔEvdw), electrostatic interaction energy (ΔEele), and internal energy (ΔEint). ΔGsol is the combination of ΔGpolar and ΔGnonpolar, which represent the polar and nonpolar solvation contributions, respectively.where the ΔGpolar term was calculated using the Poisson–Boltzmann (PB) model with the solute and solvent dielectric constants of 1.0 and 80.0, respectively. Moreover, the grid space was set to 0.5 Å. The ΔGnonpolar term was calculated by the following equation:where SASA denotes the solvent-accessible surface area, which is computed via the linear combination of the pairwise overlap approach. γ and β represent the surface tension and regression deviation of linear relations, which were set to 0.00542 kcal mol–1 Å–2 and 0.92 kcal mol–1, respectively. In addition, the residue-based energy decomposition of total energy was used to evaluate the energy contribution of each residue to the binding of RBD to the initial-type inhibitor LCB3 and mutant-type inhibitors.

Results and Discussion

Equilibrium of the Systems

The LCB3/RBD, LCB3F30Y/RBD, LCB3H6Y/RBD, LCB3L17F/RBD, LCB3H6Y-L17F/RBD, and LCB3H6Y-M7L-L17F/RBD complexes were subjected to 2 × 400 ns conventional MD simulations to adequately sample the stable complex conformational space. The root-mean-square deviation (RMSD) of the atomic position, relative to the initial structure was evaluated by least-squares fitting for the Cα atoms to investigate the overall conformational stability of the protein–ligand complex in the MD simulation. Generally, when the system obtains a stable RMSD curve, which means that it reaches an equilibrated state; the other side, the sharp fluctuations of RMSD curve indicate low structural stability. Figure , as well as Figures S1 and S2 in the Supporting Information (replica simulation) depict the smoothed fluctuation curves of the average RMSD and the time-dependent RMSD curves of the Cα atoms of the original and five mutant systems, respectively.

Figure 1

Fluctuation plot of the smooth RMSD average value of the initial complex LCB3/RBD (black) and the mutant complexes LCB3F30Y/RBD (red), LCB3H6Y/RBD (navy), LCB3L17F/RBD (dark cyan), LCB3H6Y-L17F/RBD (magenta), and LCB3H6Y-M7L-L17F/RBD (violet) along the MD simulation time, respectively. The transparent curves are the RMSD value fluctuation and the detailed RMSD fluctuation plots are shown in Figure S1 in the Supporting Information. The average RMSD curves were plotted by the Savitzky–Golay method with the adjacent 1000 data points.

As illustrated in Figure , as well as Figures S1A–C, for LCB3/RBD and mutant F30Y and H6Y complexes, the RMSD curves show that the overall structure of the complexes is in a stable state throughout the simulation process, with an average RMSD value of 1.4 Å. In addition, the mutants L17F and H6Y-L17F maintain an average RMSD value of 1.5 Å during the simulation process after 12 and 80 ns, respectively (Figure , as well as Figures S1D and S1E). However, compared with other complexes, the RMSD value of mutant H6Y-M7L-L17F shows a large fluctuation at 175 ns in the MD simulation. This is caused by the change of loop and part of the helical bundle regions in the original binding interface (see Figure S3 in the Supporting Information). After that, the RMSD value remains in equilibrium until the end of the MD simulation, with an average of 1.7 Å (Figure and Figure S1F).

Moreover, the clustering analysis was also performed on the basis of the Cα atoms of each RBD-inhibitor system to provide evidence for the structural stability of the protein. Clustering analysis was performed based on MD trajectory after 10 ns, because of the large structural changes in the early stages of the trajectory. In the large clusters composed of all trajectories at production stages, the average RMSD distance between points in each cluster is <1.7 Å, which shows that the conformational sampling is relatively concentrated in the simulation process. In this case, it has been provided that the conformation of the system is in a dynamic equilibrium state under the simulation conditions. Therefore, the conventional MD simulation trajectories are suitable for further energy analysis.

Structural Flexibility of the Systems

To investigate the structural flexibility in the entire simulation, the root-mean-square fluctuation (RMSF) of the main chain Cα atoms, reflecting the fluctuation of the atoms from their average position, was calculated on a stable trajectory. Generally, the higher value of the RMSF implies that the residue has larger flexibility. Similarly, the lower value of the RMSF means the restriction of residue movement, thereby reducing their flexibility.

As shown in Figure , for the LCB3/RBD system, residues 366–370 in the helix region and residues 371–373, 445–446, 476–481, and 518–520 in the loop region of RBD show large fluctuations with RMSF values larger than 1.5 Å, indicating that these regions are flexible. The F30Y, H6Y, H6Y-L17F (Figure C), and H6Y-M7L-L17F mutations in the five mutants result in small fluctuations of residues 444–449 in the loop of RBD relative to the initial system, and its RMSF value decreases from 1.8 Å to 1.0–1.2 Å, whereas the L17F mutation results in larger fluctuations in the above region, relative to the initial system, with the RMSF value increasing from 1.8 Å to 2.3 Å (Figure B). The H6Y mutation also causes greater fluctuation in residue 468, compared to other systems. Moreover, the L17F mutation results in a higher fluctuation of residues 369–373 in the helix and loop regions of RBD, relative to the LCB3/RBD system. However, the fluctuation in the H6Y-L17F mutation is relatively small in the above-mentioned region. Meanwhile, the L17F and H6Y-L17F mutations also exhibit synchronously larger fluctuations in the residue range 476–481 of the RBD relative to the LCB3/RBD system (Figures B and 2C). Note that the RMSF values of all residues of all inhibitors show small fluctuations, and their values are all <2.0 Å, which is consistent with the RMSD analysis, indicating that the inhibitor ligands are in a stable state and therefore exhibit lower RMSF values.

Figure 2

(A) Plots of root-mean-square fluctuation (RMSF) of the Cα atoms of the mutant complexes: LCB3F30Y/RBD (orange), LCB3H6Y/RBD (violet), LCB3L17F/RBD (green), LCB3H6Y-L17F/RBD (magenta), and LCB3H6Y-M7L-L17F/RBD (blue), relative to the LCB3/RBD complex (black) from conventional MD simulations, respectively. The regions with large change of the RMSF values are specially marked in different colors: Residue 18–25 (violet), Residue 40–48 (red), Residue 366–373 (magenta), Residue 444–449 (blue), Residue 468 (olive), Residue 476–481 (orange), and Residue 518–520 (wine). X-axis: inhibitor residue IDs from 1 to 64 are shown in light blue, and RBD residue IDs from 334 to 527 are shown in green. (B, C) Superimposed 3D-structures of different residue regions of LCB3L17F/RBD complex (highlighted in red color) RBD or LCB3H6Y-L17F/RBD complex (highlighted in magenta color) RBD with LCB3/RBD complex (highlighted in blue color) RBD having high or low RMSF values.

Interaction analysis was performed to determine the ligand–receptor interaction and the binding strength of the inhibitors to the active sites of RBD, and then to evaluate the effect of key residues on the stability and binding of the complexes in conventional MD simulations. The hydrogen bond and salt bridge interactions formed between the ligands and the RBD and their lifetime were plotted and calculated, as shown in Table , as well as Table S2 in the Supporting Information (replica simulation) and Figure .

Table 1

Key Interactions between the RBD of SARS-CoV-2 and the Inhibitors in the Protein Complex during MD Simulationsa

		Occupancy (%)
donor	acceptor	LCB3/RBD	LCB3F30Y/RBD	LCB3H6Y/RBD	LCB3L17F/RBD	LCB3H6Y-L17F/RBD	LCB3H6Y-M7L-L17F/RBD	interaction
Asn1	Thr500	53.90	50.80	48.05	55.25	43.93	0.27	H-bond
Asn1	Asn501	52.15	58.45	53.55	71.66	44.90	8.57	H-bond
Tyr14	Asp420	86.73	88.81	72.79	85.31	71.74	83.08	H-bond
Gly502	Glu4	89.73	86.53	84.53	88.93	90.23	34.01	H-bond
Arg403	Asp11	22.16	20.96	60.72	213.47	70.04	137.61	salt bridge
Tyr449	Asp3	87.78	124.49	129.49	109.45	2.65	72.84	H-bond
Lys417	Asp11	67.29	67.04	69.09	69.74	61.04	68.24	salt bridge
Tyr489	Glu34	95.33	93.80	77.29	67.59	57.42	56.10	H-bond
Gln498	Asp3	68.87	29.71	38.88	52.27	101.80	8.57	H-bond
Asn487	Glu34	7.52	20.39	2.10	0.02	1.15	31.21	H-bond
Tyr6	Gln493	–	–	40.05	–	108.00	6.30	H-bond

The hydrogen bonds and salt bridges were determined by the distance between the donor and acceptor atoms must be ≤3.5 Å and the supplementary angle of the angle formed by the donor, hydrogen, and receptor atoms must be ≤30°. Occupancy is the percentage of time that the hydrogen bond and salt bridges were formed during the MD simulation time.

Figure 3

Detailed interactions between the RBD of SARS-CoV-2 and the inhibitors (A) LCB3, (B) LCB3F30Y, (C) LCB3H6Y, (D) LCB3L17F, (E) LCB3H6Y-L17F, and (F) LCB3H6Y-M7L-L17F. The different regions of the receptor and the ligand are represented by different colors. Residues with the same color represent that they are located in the same region whereas residues with different colors represent that they are located in different regions. The key residues that form hydrogen bonds and salt bridges are shown with the stick model. Hydrogen bonds and salt bridges are represented with green dashed lines. For all interaction types, those with a probability of <30% in the trajectories of simulation are not marked in the figure.

For the complex systems including LCB3/RBD, LCB3H6Y/RBD, LCB3F30Y/RBD, LCB3L17F/RBD, and LCB3H6Y-L17F/RBD (Figures A–E), it can be found that the residue Asn1 at the N-terminal of the inhibitors forms a hydrogen bond interaction with the side chains of the residues Thr500 and Asn501 on the other side, respectively. However, there is no such interaction in the LCB3H6Y-M7L-L17F/RBD complex (Figure F). The residue Tyr14 on the first helix bundle of the inhibitors contacts with the carbonyl oxygen of Asp420 by a hydrogen bond, which exists in all complex systems. The H-bond has the highest occupancy in the LCB3F30Y/RBD complex system with an occupancy percentage of 88.81%. In all systems, the side chain of Glu34, the only residue on the second helix bundle that interacts with RBD, forms a hydrogen bond with Tyr489. In addition, Tyr449 hydrogen-bonds with the carbonyl oxygen on the side chain of Asp3, with an occupancy of 129.49%, 124.49%, and 109.45%, which means two hydrogen bonds are formed simultaneously in the same trajectory frame in the complexes LCB3H6Y/RBD, LCB3F30Y/RBD, and LCB3L17F/RBD, respectively.

Note that, compared with other inhibitor–RBD complex systems, the hydrogen bond interaction formed between Gln498 and Asp3 exhibits the highest occupancy percentage of 101.80% in the complex LCB3H6Y-L17F/RBD. Furthermore, because of the mutation of histidine residue into tyrosine residue, a new hydrogen bond is established between the side chains of Tyr6 and Gln493, which is only present simultaneously in the LCB3H6Y/RBD complex. At the same time, residue Gly502 interacts with the side chain of Glu4 via a hydrogen bond, showing the highest occupancy percentage of 90.23% in the LCB3H6Y-L17F/RBD complex. All of these indicate that there is a reason for LCB3H6Y-L17F to be an inhibitor with higher inhibitory potential than the initial inhibitor LCB3 and single-point mutant inhibitors.

Remarkably, the guanidine group of residue Arg403 on the RBD forms a salt bridge with the side chain of Asp11. This interaction is rarely formed in LCB3/RBD and LCB3F30Y/RBD systems, but it exhibits a strong interaction in LCB3L17F/RBD and LCB3H6Y-M7L-L17F/RBD complexes. The occupancies of up to 213.47% and 137.61% mean that this interaction is essential for binding and inhibiting the RBD. Although positively charged Lys417 interacts with negatively charged Asp11 in all complex systems by forming a salt bridge, it shows a relatively high occupancy of 69.74%, 69.09%, and 68.24% in the LCB3L17F/RBD, LCB3H6Y/RBD, and LCB3H6Y-M7L-L17F/RBD complexes, respectively. In addition, the hydrogen bond interaction formed between Asn487 and Glu34 shows the optimal occupancy in the LCB3H6Y-M7L-L17F/RBD system, relative to the other systems. The above also provides a basis for LCB3H6Y-M7L-L17F as the “best” inhibitor, which corresponds to the results of the binding free-energy calculation (see Table , presented later in this work).

Table 3

Binding Free-Energy Results of the Inhibitor–RBD Complexes Obtained Using MM/PBSA Calculationa

system	LCB3/RBD	LCB3F30Y/RBD	LCB3H6Y/RBD	LCB3L17F/RBD	LCB3H6Y-L17F/RBD	LCB3H6Y-M7L-L17F/RBD
van der Waals energy	–88.5 ± 0.1	–86.1 ± 0.1	–85.7 ± 0.1	–91.4 ± 0.1	–90.8 ± 0.1	–70.7 ± 0.2
electrostatic energy	–192.3 ± 0.2	–195.2 ± 0.2	–187.4 ± 0.2	–204.6 ± 0.3	–197.2 ± 0.2	–189.9 ± 0.3
polar solvation energy	221.4 ± 0.4	220.9 ± 0.4	209.3 ± 0.4	232.1 ± 0.4	223.5 ± 0.3	193.4 ± 0.4
SASA energy	–10.2 ± 0.0	–9.9 ± 0.0	–10.0 ± 0.0	–10.6 ± 0.0	–10.4 ± 0.0	–8.9 ± 0.0
binding free energy	–69.7 ± 0.3	–70.3 ± 0.2	–73.8 ± 0.3	–74.4 ± 0.3	–74.9 ± 0.3	–76.1 ± 0.3
ΔΔGbind	0	–0.6	–4.1	–4.7	–5.2	–6.4

The relative binding free energy: ΔΔGbind = ΔGmutant – ΔGLCB3/RBD. All energies are given in units of kcal/mol.

Since hydrogen bonds are one of the essential intermolecular interactions and play a crucial role in the binding of inhibitors and RBD protein, the number of hydrogen bonds between inhibitors and RBD protein was also counted, as shown in Figure S4 in the Supporting Information. The average number of hydrogen bonds between the original inhibitor LCB3 and the mutant inhibitors LCB3F30Y, LCB3H6Y, and LCB3H6Y-L17F and RBD are very similar, while the number of hydrogen bonds between the mutants LCB3L17F and LCB3H6Y-M7L-L17F and RBD is more variable throughout the simulations.

To elucidate whether these designed inhibitors could be adapted to the RBD protein in the binding interface, the SARS-CoV-2 RBD-ACE2 complex (PDB ID 6m0j(75)) was also performed with 200 ns MD simulation to determine the interface residues between RBD and ACE2. The RMSD analysis indicates that the structure has reached an equilibrium state (see Figure S5 in the Supporting Information). As shown in Table , the interactions between RBD and ACE2 are dominated by hydrogen bonds and salt bridges. Among them, the donor residues Tyr505, Gly502, Thr500, Gln498, and Gln493 on RBD interact with residues Glu37, Lys353, Asp355, Asp38, and Glu35 on ACE2 through hydrogen bonds, respectively. In contrast, residues Gln498 and Gln496, Gln493, Asn487, and Ala475 on the RBD act as receptors to form hydrogen bond interactions with residues Lys353, Lys31, Tyr83, and Ser19 on ACE2, respectively. In addition, Lys417 on the RBD forms a salt bridge interaction with Asp30 on ACE2. A comparison of these above interactions with RBD-inhibitor interactions (Table ) reveals that residues on these inhibitors can not only form hydrogen bonds and salt bridges with residues on the RBD, interacting with ACE2 residues, such as Gly502, Thr500, Gln498, Gln493, Asn487, and Lys417, but also interact with other residues on the RBD, such as Asn501, Tyr489, Tyr449, Asp420, and Arg403 through hydrogen bonds or salt bridges. Therefore, the results of the comparison suggest that these inhibitors have adapted to the RBD protein in the binding interface, further forming stable complexes.

Table 2

Key Interactions between the RBD of SARS-CoV-2 and the ACE2 in the Protein Complex during MD Simulationsa

donor	acceptor	occupancy (%)	interaction
Tyr505	Glu37	92.95	H-bond
Gly502	Lys353	92.16	H-bond
Thr500	Asp355	82.92	H-bond
Gln498	Asp38	95.76	H-bond
Gln493	Glu35	111.58	H-bond
Lys417	Asp30	94.01	salt bridge
Lys353	Gln498	55.57	H-bond
Lys353	Gln496	27.09	H-bond
Tyr83	Asn487	82.59	H-bond
Lys31	Gln493	45.24	H-bond
Ser19	Ala475	52.41	H-bond

Occupancies of >20% are shown in the table.

Dynamic Correlations between Inhibitors and the RBD

In order to further investigate the conformational changes and the correlation in the movement of protein residues after inhibitor mutation, dynamic cross-correlation matrix (DCCM) analysis was performed on the basis of Cα atoms throughout the MD simulation (see Figure and Figure S6 in the Supporting Information (replica simulation)). In the figure, the dense color regions represent the respective positive and negative correlation motion between the protein residues, in which the correlation motions (highly positive) of residues are represented as cyan color regions, whereas the anticorrelation motions (highly negative) of residues are represented as magenta color regions.[76] The uncolored regions indicate no correlation in the movement of protein residues. The highly correlated motions mean the strong interaction of the inhibitor and its mutants with the RBD binding sites, resulting in coordinated movement of the entire structure.

Figure 4

Dynamic cross-correlation map of the Cα atoms of the (A) LCB3/RBD and mutants (B) LCB3F30Y/RBD, (C) LCB3H6Y/RBD, (D) LCB3L17F/RBD, (E) LCB3H6Y-L17F/RBD, and (F) LCB3H6Y-M7L-L17F/RBD complexes. X-axis: inhibitor residue IDs from 1 to 64 are shown in light blue, and RBD residue IDs from 334 to 527 are shown in green. Y-axis: key residue IDs are additionally labeled in panel (F) for the LCB3H6Y-M7L-L17F/RBD complex.

As shown in Figures A–C and 4F, as a whole, the cyan and magenta color regions of LCB3F30Y/RBD, LCB3H6Y/RBD, and LCB3H6Y-M7L-L17F/RBD are larger and more intense than those of LCB3/RBD, which means that the correlation or anticorrelation movement of LCB3F30Y, LCB3H6Y, and LCB3H6Y-M7L-L17F is enhanced by mutation. In addition, the region marked with a black circle shows a clear difference between the mutant systems and the initial system. In LCB3H6Y-M7L-L17F/RBD, part of the 3-helix inhibitor shows a significant negative correlation relative to the inhibitor itself, which is consistent with the large fluctuations between the residues 18–25 and 40–48 of the inhibitor in the RMSF analysis. At the same time, a strong negative correlation movement is also reflected in this region of the LCB3F30Y/RBD, LCB3H6Y/RBD, and LCB3H6Y-L17F/RBD systems (Figures B, 4C, and 4E). The above analysis shows that the occurrence of residue mutations on the inhibitor LCB3 can significantly change the local dynamic characteristics, especially in the inhibitor 3-helix region.

Note that the region marked with a black box shows the correlation movement between the inhibitor and RBD. In particular, in the LCB3H6Y-M7L-L17F/RBD (Figure F), the inhibitor residues 10–30 exhibit strong negative correlation movement with residues 401–434 and 460–468 of the RBD. Also, the negative correlations between residues 48–64 and 30–48 of the inhibitor and residues 401–434 and 434–455 of the RBD, respectively, increased to some extent, relative to the other systems. At the same time, strongly negative correlation movements are also reflected between residues 30–58 of the inhibitor and residues 492–511 of the RBD. Thus, the DCCM analysis of the inhibitor LCB3H6Y-M7L-L17F bind RBD residues reveals the strongest correlation among the residues movements, which is attributed to the most robust interaction.

Dominant Motions of Inhibitor Ligands

The dominant motions of inhibitor ligands were studied to evaluate their dynamic coordination with RBD. The porcupine plots of the first eigenvector (PC1) for the ligands and RBD were generated to qualitatively monitor the difference in the dominant motions of the five mutant inhibitor ligands and the initial inhibitor ligand.[77] As shown in Figure , the length of the arrow indicates the amplitude of the movements of the complexes, and the cone indicates the direction of the movement of the complexes. The mutants LCB3H6Y-L17F and LCB3H6Y-M7L-L17F present strong motion relative to the other ligands. Especially in the LCB3H6Y-M7L-L17F system, because the inhibitor exhibits a movement with the backbone pointing to the RBD direction. This bottleneck points to the binding site of the inhibitor and RBD, thereby further promoting the binding of the inhibitor and RBD. The single-point mutants LCB3F30Y, LCB3H6Y, and LCB3L17F also present different motion intensity and amplitude from LCB3. These collective motions are due to the strong and stable interaction between the inhibitor ligands and the RBD. Thus, from this perspective, it can be found that the specific interactions between miniproteins and RBD (Table ) can cause some slight differences in motion, whereas the combination of all these interactions results in the synchronous motions of the receptor–ligand as a whole.

Figure 5

Representation of the dominant motions of the residues obtained from the first principal component (PC1) for the complexes established between RBD (sky blue) and (A) LCB3 (cyan), (B) LCB3F30Y (lime green), (C) LCB3H6Y (salmon), (D) LCB3L17F (pink), (E) LCB3H6Y-L17F (orange), and (F) LCB3H6Y-M7L-L17F (pale green).

Binding Free-Energy Calculation

The Molecular Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) method is a common approach used to calculate the binding free energy of small ligands to biological macromolecules.[78] The method cannot accurately reproduce the absolute experimental binding free energy, but it can give a good rank,[79] and it is used to evaluate the binding ability of inhibitor molecules to RBD. In this study, this method was implemented to calculate the total binding free energy of RBD with initial and mutant inhibitors. In order to clarify the qualitative effect of mutations on the binding affinity of RBD, a total of 6000 structures from the MD trajectory with stable conformations were selected to calculate the contribution of individual components such as ΔEvdw, ΔEele, ΔEpolar, and ΔEnonpolar, to obtain the total binding free energy ΔGbind. The prediction results of binding free energy and each energy term are shown in Table , as well as Table S5 in the Supporting Information (replica simulation). Thereinto, the negative value means that the interaction is favorable whereas the positive value indicates that the interaction is unfavorable. According to the simulation results, compared with the interaction between LCB3 and RBD (−69.7 ± 0.3 kcal/mol), the interaction of inhibitors LCB3F30Y, LCB3H6Y, LCB3L17F, LCB3H6Y-L17F, and LCB3H6Y-M7L-L17F compounds with RBD is favorable with ΔGbind being −70.3 ± 0.2 kcal/mol, −73.8 ± 0.3 kcal/mol, −74.4 ± 0.3 kcal/mol, −74.9 ± 0.3 kcal/mol, and −76.1 ± 0.3 kcal/mol, respectively.

In the LCB3/RBD system, the triple helix ligand as an inhibitor of RBD is used to provide a reference for the design of new inhibitors in this study. On this basis, compared with the LCB3/RBD complex system, the LCB3L17F/RBD system in all the single-point mutants possesses the largest relative binding free energy, −4.7 kcal/mol, which means that the inhibitory constant (K) of the inhibitor LCB3L17F decreases to ∼1/2655 of the LCB3 inhibitory constant (ΔG = RT ln K), and the potency of inhibiting RBD is significantly increased in comparison with LCB3. This is consistent with the enhanced receptor–ligand interaction in the interaction analysis. Second, for the inhibitors of double-point mutation LCB3H6Y-L17F, the calculated value of the relative binding energy (ΔΔGbind = −5.2 kcal/mol) indicates the inhibitory activity on RBD is further enhanced. The “best” inhibitor in all systems, according to the relative binding energy calculation results, is LCB3H6Y-M7L-L17F, which has the lowest relative binding energy value of −6.4 kcal/mol. Therefore, the five mutant inhibitors—LCB3F30Y, LCB3H6Y, LCB3L17F, LCB3H6Y-L17F, and LCB3H6Y-M7L-L17F—are all potential inhibitors of RBD, and their inhibitory abilities are successively increased.

Further analysis of the calculated energy component of the binding affinity shows that the contributions of the van der Waals interaction, electrostatic interaction, and nonpolar solvation energy work together to facilitate the formation of the complexes in all systems. As shown in Table , from these values, electrostatic interaction seems to be the main reason for the formation of inhibitor–RBD complexes, which are −192.3 ± 0.2, −195.2 ± 0.2, −187.4 ± 0.2, −204.6 ± 0.3, −197.2 ± 0.2, and −189.9 ± 0.3 kcal/mol, respectively. The contributions of the calculated nonpolar solvation energy based on the solvent accessible surface area (SASA) are also favorable for the formation of the complexes, although their contribution is small. Furthermore, the triple mutated LCB3H6Y-M7L-L17F/RBD complex can exhibit the lowest binding free energy, which is critically related to the polar solvent energy with relatively less unfavorable contribution, compared to other systems.

Per-Residue Energy Decomposition Analysis

The residue-based energy decomposition of ΔGbind is an effective analysis method to evaluate the energy contribution of key residues. The energy contributions of these residues to the total binding energy are shown in Table . The critical residues with favorable energy contribution mainly include the acidic amino acids of inhibitors, such as aspartic acid and glutamic acid, and the basic amino acids of RBD, such as arginine and lysine. These residues are generally involved in the electrostatic interaction during the MD simulation process, which means that the electrostatic interaction drives the binding. Then, the difference in energy contributions of each residue between the mutant and initial systems (ΔΔG = ΔGmutant – ΔGLCB3/RBD) was plotted to clarify the key residues, as shown in Figure .

Table 4

Residue-Based Energy Decomposition of Key Favorable Contributing Residues and Residues at Mutation Sites for Inhibitors LCB3, LCB3F30Y, LCB3H6Y, LCB3L17F, LCB3H6Y-L17F, and LCB3H6Y-M7L-L17F and RBDa

residue	LCB3/RBD	LCB3F30Y/RBD	LCB3H6Y/RBD	LCB3L17F/RBD	LCB3H6Y-L17F/RBD	LCB3H6Y-M7L-L17F/RBD
D2	–12.0	–11.8	–11.4	–11.3	–12.7	–11.8
D3	–8.8	–10.0	–10.0	–9.5	–10.4	–11.8
E4	–13.2	–12.9	–13.0	–12.9	–13.0	–14.4
H/Y6	1.8	1.2	0.1	1.6	0.9	–0.2
M/L7	–4.0	–3.9	–3.5	–4.0	–3.7	–2.4
L8	–1.3	–1.3	–1.2	–1.4	–1.3	–0.7
D11	–9.9	–10.5	–8.1	–8.4	–7.6	–6.2
T14	–2.3	–2.3	–2.2	–2.2	–2.6	–1.9
E15	–11.5	–11.6	–11.7	–11.9	–11.7	–11.6
L/F17	–0.6	–0.5	–0.6	–1.1	–1.2	–0.8
D22	–6.0	–6.1	–5.8	–6.1	–6.1	–5.9
E23	–6.0	–7.0	–5.7	–6.1	–5.9	–4.9
E24	–5.7	–5.7	–5.5	–5.7	–5.7	–5.5
F/Y30	–2.1	0.2	–2.0	–2.3	–2.2	–2.0
F33	–1.6	–1.5	–1.5	–1.6	–1.8	–1.7
E34	–3.9	–3.1	–4.0	–4.4	–4.6	–2.5
D37	–7.2	–6.6	–6.4	–6.5	–7.9	–6.6
T40	0.2	–0.9	–1.0	–0.6	–1.1	–0.4
D44	–7.3	–7.0	–6.9	–6.9	–8.0	–6.9
E49	–8.9	–8.9	–8.9	–8.8	–9.3	–8.7
E53	–7.9	–7.9	–8.0	–8.1	–8.0	–7.8
E54	–6.6	–6.6	–6.6	–6.7	–7.0	–6.6
E57	–6.4	–6.5	–6.4	–6.6	–6.7	–6.5
E60	–6.9	–7.2	–6.9	–7.2	–7.0	–7.1
S64	–5.5	–5.8	–5.3	–5.7	–6.3	–5.7
N334	–6.6	–6.6	–6.5	–6.5	–6.6	–6.5
R346	–10.6	–10.4	–10.3	–10.4	–10.7	–10.5
R355	–9.5	–9.6	–9.4	–9.5	–9.4	–9.5
K356	–8.8	–8.7	–8.6	–8.8	–8.6	–8.7
R357	–7.7	–7.8	–7.6	–7.7	–7.5	–7.6
K378	–13.5	–13.6	–13.1	–13.2	–13.3	–13.3
K386	–7.7	–7.8	–7.6	–7.6	–7.7	–7.7
R403	–11.1	–11.8	–11.8	–10.0	–11.4	–13.2
R408	–21.8	–22.5	–21.8	–21.8	–21.9	–22.5
K417	–6.3	–6.9	–7.1	–6.9	–7.4	–8.7
K424	–12.4	–13.0	–12.4	–12.5	–12.5	–12.6
K444	–12.7	–12.8	–12.3	–12.3	–13.1	–12.4
R454	–11.4	–11.8	–11.7	–11.7	–11.6	–11.8
L455	–2.1	–1.9	–2.0	–2.0	–2.0	–1.6
F456	–1.7	–2.1	–1.6	–1.7	–1.7	–1.5
R457	–10.7	–11.1	–10.8	–10.6	–10.4	–10.8
K458	–10.4	–10.8	–10.3	–10.6	–10.2	–10.1
K462	–9.2	–9.4	–9.0	–9.1	–8.9	–9.0
R466	–10.2	–10.4	–10.1	–10.3	–10.0	–10.2
Y505	–1.6	–1.3	–1.8	–1.7	–1.4	–0.5
R509	–11.3	–11.4	–11.2	–11.4	–11.1	–11.6

Energy values are given in units of kcal/mol and energy values above −1.0 kcal/mol are considered. See Table S6 in the Supporting Information for residue-based energy decomposition on most of the residues in all systems. Inhibitor residue IDs from 1 to 64 are shown in the upper part of the table, and RBD residue IDs from 334 to 527 are shown in the lower part of the table.

Figure 6

Energetic differences of the residue contributions to the binding free energies between the initial system and mutant system (ΔΔG = ΔGmutant – ΔGLCB3/RBD). (A) LCB3F30Y + RBD and LCB3 + RBD (cyan), (B) LCB3H6Y + RBD and LCB3 + RBD (green), (C) LCB3L17F + RBD and LCB3 + RBD (blue), (D) LCB3H6Y-L17F + RBD and LCB3 + RBD (olive), and (E) LCB3H6Y-M7L-L17F + RBD and LCB3 + RBD (violet). X-axis: inhibitor residue IDs from 1 to 64 are shown in light blue, and RBD residue IDs from 334 to 527 are shown in green.

The mutation of phenylalanine to tyrosine in the F30Y mutant results in the conversion of its own energy contribution to a disadvantage. However, the energy of other surrounding residues Asp3, Asp11, Glu23, and RBD residues Arg403, Arg408, Lys417, and Lys424 have significantly increased (Figure A). The above-mentioned residues are either acidic amino acids or basic amino acids, which is why electrostatic interaction is the main driving force for binding. Simultaneously, the negative contributions of His6, Lys41, and Ser494 to the total binding energy were weakened (Table S6), whereas Tyr40 has changed from unfavorable to favorable. Furthermore, the energy contribution of the individual residues of Glu34 and Asp37 is slightly decreased, but they still promote the overall binding. However, the contribution of the residues Asp405 and Asp420 on the RBD to the total binding energy is found to be more unfavorable. In summary, the mutant LCB3F30Y does not significantly increase the binding to RBD in comparison with LCB3.

In LCB3H6Y/RBD system, the residue His6 of the ligand is mutated to tyrosine. Although the contribution of the mutated Tyr6 to the total binding energy remains unfavorable, the binding energy is significantly reduced by −1.7 kcal/mol (see Table and Figure B), this situation also occurs in the inhibitor residue Lys41 and RBD residues Tyr453, Gln493, and Ser494. In addition, the mutation changes the contribution of residue Tyr40 to the total binding energy in a favorable direction. Moreover, residues that are originally favorable for binding, such as Asp3, Arg403, and Lys417, become more favorable for ΔGbind, which is consistent with the enhanced interaction between Arg403–Asp11, Lys417–Asp11, and Asp3–Tyr449 in the interaction analysis. Although the contribution of residues Asp2, Asp11, and Asp37 to the total binding energy is decreased, they still have a beneficial effect on the binding of the inhibitor to the RBD protein. The residues Glu406 and Glu484 of RBD are more unfavorable for binding. But overall, the inhibitory potential of mutant inhibitor LCB3H6Y on RBD is further improved, compared with LCB3.

The mutation of leucine to phenylalanine in the L17F mutant causes a substantial increase in the individual energy contribution of residue Phe17, which provides more favorable energy for binding than the initial system. In addition, the residues Asp3, Met7, Glu15, Asp22, Glu23, Glu24, Phe30, Glu53, Glu54, Glu57, Glu60, and Ser64 with beneficial energy contributions on the inhibitor have different degrees of energy increase (Table ). Furthermore, Figure C shows that the RBD residues Lys417, Phe486, and Asn487 have a larger energy increase in their contribution to the total energy. Meanwhile, the unfavorable contribution from inhibitor residues Asn1, Thr10, Lys27, and Lys41 is reduced to some extent, which also occurs on the residues Asp405 and Ser494 of RBD. Moreover, although the contribution of the residues Asp2, Asp11, and Asp37 on the inhibitor, and Arg403 on RBD to the total binding energy decreases slightly, they still play a critical role in binding, because of their high contribution value. In other words, the mutation makes almost all residues more conducive to binding, so the inhibitory ability of the inhibitor LCB3L17F on RBD is further improved compared to LCB3.

For the H6Y-L17F mutant, the simultaneous mutation of two residues results in a certain increase in the energy contribution of residues Tyr6 and Phe17. It is also found that compared with the initial inhibitor, most of the residues with favorable contributions, such as Asp2, Asp3, Thr14, Glu15, Asp22, Glu24, Phe30, Phe33, Glu34, Asp37, Tyr40, Asp44, Glu49, Glu53, Glu54, Glu57, Glu60, and Ser64, have increased to varying degrees (Table ), which can be explained by the enhanced interaction caused by the residue mutation. Moreover, it can be observed from Figure D that the unfavorable contribution of residue Lys27 to the total binding energy is reduced. Furthermore, the mutation also increases the energy contribution of residues Lys417, Gly446, and Phe486 on the RBD, at the same time decreasing the energy of the unfavorable contribution of residues Asn487 and Ser494, thereby turning the unfavorable contributions of residues Gly447 and Tyr449 into favorable contributions. Since the energy increase of residues favorable for binding far exceeds the unfavorable increase. Therefore, the double-point mutant inhibitor LCB3H6Y-L17F shows better inhibitory potential than the single-point inhibitor.

Based on the H6Y-L17F mutant, to further improve the ability to inhibit RBD. When three different residues are simultaneously mutated (Figure E), the residues Tyr6, Tyr40, and Asn43 all change from an unfavorable contribution to binding to a favorable contribution, and the energy contribution of residue Phe17 at the mutation point increases, while the energy contribution of Leu7 decreases slightly. At the same time, the inhibitor residues Asp3 and RBD residues Arg403 and Lys417 exhibit a significant increase in energy to facilitate binding, which is consistent with the enhancement of the salt bridge interaction between them. Figures D and 6E shows that LCB3H6Y-M7L-L17F has more individual residues with more favorable energy contributions for total energy than LCB3H6Y-L17F, the double-point mutation inhibitor. In particular, such mutations greatly increase the beneficial energy contribution of the terminal residues of the inhibitor such as Asn1, Asp3, Glu4, and Tyr6, which also makes the residues Arg403, Arg408, Lys417, and Phe486 on the RBD contribute to the total binding energy contribution further improve, and the residues Tyr453, Asn487, and Tyr495 with unfavorable contributions are converted into residues with favorable contributions. Meanwhile, the unfavorable contribution energy to the total binding energy of residues Gln493, Ser494, and Asn501 is reduced. However, although the contribution of residues Leu7, Leu8, Asp11, Glu23, Glu34, and Asp37 on the inhibitor to overall energy is also decreased, they still have a high negative energy value and play a positive role in the binding of the inhibitor to RBD. Therefore, LCB3H6Y-M7L-L17F is the “best” inhibitor in all systems, with the highest inhibitory potential.

Potential Uncertainties from Related Experiments

So far, the current simulations have presented several better computer-designed miniproteins to target the RBD of SARS-CoV-2. From an energetic point of view, all the inhibitors demonstrate better antiviral activity and inhibitory capacity than the original peptide, which is attributed to enhanced receptor–ligand interactions. The next work is to validate the binding affinity of these in-silico designed miniproteins to the RBD by well-designed experiments. For in vitro experiments, the affinity result is expected to be consistent with the present theoretical calculation under the same conditions, because our simulations are based on the existing experiments[44] and are in good agreement with the experimental results. For in vivo experiments, especially for some clinical experiments, the therapeutic effect of the present miniprotein inhibitors theoretically predicted should be carefully examined and determined because the complexity of human environment may be very different from the solvent condition in our simulation, which may result in some unexpected changes of the interaction between the miniproteins and the RBD.

Conclusions

In this work, we designed the five mutant inhibitors, including LCB3F30Y, LCB3H6Y, LCB3L17F, LCB3H6Y-L17F, and LCB3H6Y-M7L-L17F, both single-point, double-point, and triple-point mutants based on the initial inhibitor LCB3 to against SARS-CoV-2 spike protein. MD simulations were performed to investigate the dynamic interactions between six inhibitor ligands and RBD. The simulations and calculation of MM/PBSA binding free energy accurately explain why the mutant inhibitors have a stronger inhibitory effect compared with LCB3. In particular, because of the lowest binding free energy and relative binding free energy (−6.4 kcal/mol), the LCB3H6Y-M7L-L17F becomes the “best” inhibitor relative to several other inhibitors. The interaction study at the atomic level and the decomposition of the energy contribution of individual residue identify the binding mode and key residues of each inhibitor ligand to RBD and reveal the difference in binding affinity between the six inhibitors and RBD, thereby clarifying that the single, double, and triple-point mutant inhibitors have different effects on the binding of RBD. Moreover, the dynamic correlation between each inhibitor and RBD determines the effect of mutations on the movement correlation of residues on the inhibitors. In conclusion, this study clarifies the potential mechanism of the RBD contacting with inhibitors and examines the effect of mutations on inhibitory activity. Therefore, some new mutant inhibitors are proposed that will provide a better interaction mechanism for the rational design of RBD inhibitors of SARS-CoV-2.