Synergistic Effects of Natural Compounds Toward Inhibition of SARS-CoV-2 3CL Protease.

The biggest challenge in medical management and control of the COVID-19 pandemic is the nonavailability of the treatment molecules. While vaccines and other biotherapeutic products for managing COVID-19 have reached the market, a small-molecule cure is yet to be developed. This is relevant because the cost of production, storage, and ease of distribution of a small-molecule drug are significantly more favorable than those of biologics. In this paper, we present a multicompound approach, where two drug molecules are administered concurrently to offer an effective therapy for COVID-19. The co-action of the two compounds, each derived from natural origins, has been demonstrated against the 3CL protease, already recognized as a potential drug target for inhibiting SARS-CoV-2. The pair of compounds pursued in this study are flavonoid and naphthalene scaffold. Individually, they offer ∼30 to 35% inhibition at 10 μM. Comprehensive docking and molecular dynamics simulations elucidate that these compounds exhibit excellent binding in the process, which however quickly deteriorates, and the ligand is separated from the binding site. This suggests that while the ligands initially bind with the protease, they are unable to maintain it for an extended period. However, the simulation showed that a simultaneous docked complex of both the compounds together with the protein boosts the stronger binding for a sufficient time. The enzyme assay exhibited 97 and 85% inhibition activity when both compounds were used together at 100 and 50 μM, respectively. Later, a multiconcentration assay was used to determine the coinhibitory activity, and it was observed that the compounds have ∼20 to 30% inhibition activity even at lower concentrations of 0.5 and 1 μM. Surface plasmon resonance was used to measure the binding of the compounds, and when used together, the compounds had a 10-fold greater binding affinity. Thus, the results demonstrate a synergistic mechanism between the two compounds that enhances the inhibition activity against SARS-CoV-2 3CL protease.

Introduction

COVID-19 has mobilized research groups around the world to design and develop therapeutic molecules against SARS-CoV-2 to combat COVID-19. Under the “preventive medicine” approach, several vaccines have been recently brought to the market. However, vaccine development offers its own set of challenges in terms of manufacturing cost, transportation and storage, and mass vaccination, with most vaccines requiring double doses.[1] Variants of SARS-CoV-2 have been identified in many countries, and the efficacy of the various vaccines against the emerging SARS-CoV-2 variants is a question.[2,3] Despite the occurrence of several waves of COVID-19 around the world and limited vaccine availability, the value of using inhibitor molecules to reduce the disease’s burden has risen worldwide.[4] There is a need for a small-molecule treatment to effectively deal with the COVID-19 pandemic. Repurposing of known drug molecules against SARS-CoV-2 has been attempted with limited success.[5−11]

Protein targets for SARS-CoV-2 that are essential for the survival of the virus have been identified and are being used for drug discovery.[12−14] There are seven key drug targets for SARS-CoV-2, namely, (1) spike protein, (2) envelope protein, (3) membrane protein, (4) protease, (5) nucleocapsid protein, (6) hemagglutinin esterase, and (7) helicase.[15] Three-dimensional structures of many of these targets are being continuously solved to strengthen the structure-based drug design. Protease, a known drug target for RNA viruses, has been used for designing inhibitory molecules that can further target viral replication.[16] Natural compounds have shown diverse chemical scaffolds and have been successfully used in managing various diseases including viral infections. Moreover, about 45% of blockbuster drugs have originated from either natural products or their derivatives.[17−19] Recently, the structure of 3CL-protease has been solved in complex with the novel inhibitor 5,6,7-trihydroxy-2-phenyl-4H-chromen-4-one.[20] This protein structure can be used to screen and design natural molecules or their derivative showing inhibitory properties against 3CL-protease. Among the natural molecules, biflavonoids have been reported in several studies as a potential inhibitor for the SARS-CoV-2 protease.[21−25] Its applications are also reported as a safe and effective treatment for viral infections.[26,27] The protease of SARS-CoV-2 has multiple binding sites other than the active site where different ligands can bind. Applying multiple compounds on protease could enhance the inhibitory activity. Studies have reported combinatorial and synergetic action of two and more compounds as a rational drug discovery approach.[28−34] A combination of bioflavonoid compounds with other bioactive natural compounds can thus improve the overall response.

Computational techniques have assumed a significant role in the drug discovery process by slashing the research time and mitigating the risk factor.[35] Molecular docking and molecular dynamics (MD) simulation are two in silico techniques that guide rational drug design. While molecular docking is used to find the fit of the compound in the protein’s binding site,[36] MD simulations help evaluate the motion of the system with respect to each other.[37] These techniques can be used to understand the synergistic effect of two compounds against a single protein target.

In this paper, we use computations complemented with experimental validation to find a pair of natural derivative compounds that can target the SARS-CoV-2 protease simultaneously and improve the overall inhibitory activity. An in silico virtual screening is performed against protease. This resulted in the identification of a pair of compounds showing high binding affinity with the protease protein. The inhibitory effect of these compounds was further validated using a SARS-CoV-2 3CL Protease assay kit. Individual activity of the screened compounds was recorded, and the best two compounds were considered for further study. These two compounds have flavonoid and naphthalene chemical scaffolds, respectively. Individual and combined molecular docking of these two compounds indicated binding of these molecules at different binding sites of the protein. Further, MD simulations of these docked complexes were performed to understand the synergistic characteristic. Eventually, both compounds were applied simultaneously and sequentially at 10 and 100 μM concentrations in 3CL Pro inhibition assay. The use of these compounds in combination culminated in better inhibition characteristics. Later, a dose-dependent activity assay was performed to calculate the IC50 of the combined application of these two molecules. In addition to the enzyme-based assay, surface plasmon resonance (SPR) was also performed to quantify the protein–ligand binding affinity. This information would be helpful in determining the binding kinetics of individual molecules as well as combinations of molecules to the 3CL protease of SARS-CoV-2.

Results and Discussion

Screening of Natural Compounds

Virtual screening covered a broad range of molecules, including binders for 3CL protease and active compounds against SARS-CoV-2 or similar viruses. The ZINC natural compound used for screening and its pharmacophores were matched, and the structure with the lowest root mean square deviation (RMSD) values with the overall pharmacophore was selected for docking. In addition to the pharmacophore screening of the ZINC natural molecule repository, a separate docking-based screening was conducted to screen compounds from BIOFACQUIM and Chemfaces repository due to their small population of the compounds. These screenings facilitated the selection of the top 10 compounds for further validation with the early enzyme-based assay. The molecular structure and names of these screened compounds are shown in Table .

Table 1

Compounds Screened against 3CL Protease of SARS-CoV-2 Using the In Silico Technique

Early Enzyme Activity Assay

To evaluate the antiviral efficacy of the 10 compounds identified from the in silico virtual screening, they were tested at 10 mM concentration with a 3CL protease SARS-CoV-2 assay kit. Compounds 9 and 6 showed relatively higher activity and were rated as effective hits (Figure ).

Figure 1

Percentage inhibition activity of 10 screened compounds using an enzyme-based assay technique for 3CL protease.

Best Hits Docking (Individual and Combined)

Early enzyme-based assay at a 10 μM concentration indicated that compounds 6 and 9 exhibit activity against 3CL protease. Their average inhibition to the binding of natural substrates to the protein was 31.47 and 33.04%, respectively. Rest of the compounds displayed negligible inhibition and thus were not taken into further consideration. Compounds 6 and 9 were docked blindly to the protein, where the complete protein surface was used as the docking grid in AutoDock. Figure a shows the docked complexes for both the compounds in their best binding pose. In this blind docking, both these compounds bind at the catalytic site of 3CL protease. Docking produced 100 poses for each compound, and their respective binding energies were calculated. The boxplot of these binding energies of 100 poses is shown in Figure b. Energies are converted to positive values for readability purposes. Compound 6 has an average binding energy of −7.17 kcal/mol, while compound 9 has an average binding energy of −8.32 kcal/mol. Compound 9 performed better than compound 6, both in enzyme activity assay as well as in docking energy simulations. Figure c shows the protein–ligand interaction of compounds 6 and 9 with 3CL protease protein. While compound 6 forms two hydrogen bonds with GLU166 and CYS145, compound 9 has single hydrogen bonds with GLY143 and a stacking interaction with HIS41. Both the residues, HIS41 and CYS145, are present near compounds 6 and 9 in their docked complexes. The strong interaction of the catalytic diad of 3CL protease (HIS41 and CYS145) further advocates the meaningful binding of compounds 6 and 9 with the protein. Although the binding energies of both these compounds are highly promising with several relevant protein–ligand interactions, their performance in the experiment is relatively weaker. Docking and binding experiment assays cumulatively suggest that while they do bind to the protein, the binding does not last long to deliver high efficacy. Later, the combined binding activity of compounds 6 and 9 is studied using a multiple ligand simultaneous docking method.

Figure 2

(a) Binding poses of compounds 6 (green) and 9 (pink) for the best-docked pose. (b) Boxplot showing the binding energies of the top 100 binding poses of each compound achieved from AutoDock Vina. Energies are represented in positive values or better readability. (c) Protein–ligand interaction of compounds 6 and 9 with 3CL protease in their respective best-docked poses as per their binding energy scores in AutoDock. (d) Simultaneous binding of compounds 6 and 9 to the 3CL protease; compound 9 is docked first, and then compound 6 is allowed to dock using a blind docking approach. The table shows the neighboring residues under 4 Å for each compound.

In multiple docking, compound 9 is used to dock first because it exhibited strong binding activity both in docking and in enzyme assay. It binds at the catalytic site of the protein, as it does in the case of individual docking. Once compound 9 is docked, then compound 6 is blindly docked to the available protein surface where predocked compound 9 is treated as the element of the protein molecule. Here, compound 6 could not find the position at the catalytic site as it did in individual docking (Figure ) because compound 9 already occupied the catalytic site. Hence, compound 6 docked at a distinct position away from the main catalytic site. Docking of multiple ligands to 3CL protease resulted in seven docked poses with the best binding energy of −7.29 cal/mol, as shown in Figure d. Residues that are lined up around 4 Å from each ligand are tabulated in Figure d. Compound 6 is placed at a minimum of 17 Å away from the docked position of compound 9, as shown in Figure d. It has a new set of neighboring residues, which are also listed in Figure d, and it binds at the relatively open surface and thus does not have a large number of residues in its surrounding. There are only eight residues reaching its 4 Å proximity. SER139 is one of the neighboring residues of compound 6, which is closest to catalytic diad compared to any other residue.

A simultaneous docking study suggests that both the molecules can bind to the protein at the same time and that this can bring greater conformational changes in the protein structure. This can result in a synergistic mechanism that adds the binding affinity for both the compounds. Early enzyme assay suggested that both compounds 6 and 9 bind to the protein, but their binding does not translate into high inhibition. Short-lasting binding could be one explanation for the lower activity of these compounds. The binding of both compounds together with the protein can mutually enhance the period of interaction and thereby improve the inhibition activity.

Combined Assay (50 and 100 μM)

The antiviral efficacy of the compounds 9 and 6 against 3CL protease was determined using its combination at two concentrations, and they were added to the 3CL protease sequentially. The assay was followed as suggested in the previous section. The result demonstrates the synergistic effect against 3CL protease at 50 and 100 μM as compared to the individual tested compounds 9 and 6. The percent inhibition at 50 and 100 μM is 85 and 97%, respectively, as shown in Figure . The findings revealed that these molecules work well together and can have a significant effect, suggesting that their combination may be useful in the treatment of SARS-CoV-2.

Figure 3

(a) Activity of compounds 9 and 6 represent strong efficacy against the coronavirus 3CL proteases from SARS-CoV-2. The combination of compounds 9 and 6 was used at 100 and 50 μM concentrations. (b) GC376 at 100 μM was used as an internal control. Data are shown as mean ± S.D. for three technical replicates.

MD Simulation

MD simulation is a popular approach for elucidation of the interactions between the protein and the ligand molecules. A 100 ns long dynamic simulation was performed to capture the fluctuating behavior of the protein–ligand complex. Two complexes of protein and ligand considered for the simulation are (1) protein-9 and (2) protein-9-6. The best-docked pose generated during Autodock was selected for the MD simulation. Before the production phase, the entire system was soaked in a water solvent and equilibrated at a given temperature and pressure. This process was played under NVT and NPT ensembles for 1000 ps each where protein–ligand atoms are constrained to attain the equilibrium. Figure S1a,b show the potential energies of the complex’s protein-9 and protein-9-6 during the steepest descent energy minimization. Both the complexes are sufficiently minimized, and this indicates the minimum atomic steric clashes within the system. Similarly, temperature and pressure variations are demonstrated in Figure S1c–f. Temperature was set at 300 K during NVT equilibrium and showed a stable pattern around this value. The “running average” for the temperature is shown in the ‘red’ line, and it can be observed that the fluctuation around the mean is minimum. During NPT equilibration, reference pressure was set at 1 bar and both the complexes reacted comparably in the NVT ensemble. Here, the running average was more stable in the protein-9 complex. However, the pressure fluctuation pattern was under the acceptable range for both complexes. Figure S1 collectively shows that the selected docked complexes are sufficiently equilibrated under NVT and NPT ensembles, and the surrounding water medium is stabilized.

RMSD Fluctuation

After the completion of the equilibrium stage under the NVT and NPT ensembles for 1 ns each, a sufficiently long 100 ns production phase was performed. Here, no restraint was applied to the protein and ligand atoms. The “Protein-9” docked complex was simulated first, and the RMSD deviation of the compound 9 was recorded. RMSD of the ligand was calculated by fitting the protein and ligand molecule separately. Protein Cα RMSD was also calculated to determine the deviation inside the protein molecules. RMSD of Cα atoms and compound 9 for ligand molecule-based alignment was very marginal, as shown by “green” and “black” lines in Figure a.

Figure 4

(a) Root mean square deviation (RMSD) of compound 9 complexed with protease protein. Compound RMSD was calculated using ligand and protein both as reference, respectively. The RMSD unit is angstrom. (b) RMSD of compounds 9 and 6 complexed with protease protein. Compound RMSD is calculated using ligand itself as reference and protein as reference. The RMSD unit is angstrom. (c) Distance between the center of the mass of the active site and compounds 9 and 6, respectively. (d) Distance of compounds 9 and 6 with the catalytic dyad of 3CL protease (His41 and Cys145).

However, when the protein molecule was considered as the reference for alignment and the RMSD of compound 9 was calculated, an anomalous behavior was observed (red line in Figure a). Prior to 7 ns of simulation time, RMSD fluctuated under 3 Å, but then the ligand showed continuous elevation in RMSD. It reached 150 Å after 20 ns of simulation time, indicating that compound 9 has moved far away from the binding site and eventually from the protein itself. Later, the simulation trajectory was inspected using the VMD visualization tool. Visual inspection also showed that after 7 ns, the binding started losing its contacts and the molecule came out of the pocket at 10 ns. Later, this molecule disconnected completely from the protein, as shown in Figure S2. This confirmed that binding of compound 9 with the protease protein is a short-term binding, and it does not persist for long, as supported by the low inhibition observed in the enzyme-based assay.

Similarly, a combined docked complex of protease with both compounds 9 and 6 was studied for the RMSD for the bound molecules. Cα RMSD of the protein was calculated and shown as the “yellow” line in Figure b. Deviations in protein atoms were marginal during the simulation and stabilized at 2.5–3.0 Å from the initial conformation. Ligand’s RMSD was calculated using protein and ligand independently as reference fitting molecules. When the ligand molecule was taken as a reference, then neither compound 9 nor compound 6 produced high structural deviation from their respective starting conformation, as shown by “black” and “green” lines in Figure S2, respectively.

The magnitude of these RMSD values matches with the protein Cα RMSD. In contrast, when the protein molecule was taken for superposition and RMSD of ligands was recorded, the values were relatively higher. Here, compound 6 reached 10 Å in a short time and stabilized there throughout the simulation (blue line). This shows the initial movement of compound 6 in the binding site, followed by acquisition of a stable geometrical and thermodynamical state. RMSD of compound 9 showed a different behavior compared to that of compound 6. It stabilized to 6 Å during the first 70 ns of the simulation and then it reached 9 Å for the rest of the simulation (red line). Both the compounds changed their conformation from the initial state but stabilized with the new conformation. Compound 9 has a two-phase stability pattern where the first phase corresponds to the first 70 ns of simulation time, while the second phase covers the last 30 ns. However, compound 6 has a single-phase stability pattern that continues from 10 ns till the end of the simulation.

Distance Matrix

The shift in geometry and location of both ligands were measured relative to the active site of a protein and catalytic dyad. The distance between the centers of the mass of these ligands with the active site is shown in Figure c. Compound 9 is sitting within the active site, and its distance showed a consistent nature of 2 Å during the complete simulation, while compound 6 is docked away from the core active site and its distance relatively fluctuates higher in the range of 7–10 Å. The distance between these two ligands was also calculated (black), as shown in Figure c. His41 and Cys145 are known as the catalytic dyad for 3CL protease [28]. The distance between the center of mass of the catalytic dyad with the docked compounds is shown in Figure d. Here, compound 9 again showed a consistent and lower distance with the catalytic dyad at 2 Å, while compound 6 exhibited higher fluctuation and magnitude at 15 Å. Figure c,d shows that compound 6 is constantly placed away from the active site and catalytic dyad, and the minimum distance it showed in certain frames is ∼5 Å. However, compound 9 is always positioned within the active site of the protein and near the catalytic dyad.

Clustering

Clustering was performed under 3 Å RMSD criteria on the complete complex for the 100 ns trajectory to find the most dominant docked conformation. This resulted in the formation of 4 clusters where the top cluster had 9965 structures and the other three clusters had 21, 8, and 7 structures. This suggests that simulation produced similar conformations and the deviations are under the acceptable range.

Binding Energy

The MMPBSA method was used to calculate the binding free energy between protein and ligand in their bound state. “g_mmpbsa” is a tool developed to calculate the MMPBSA energy for the simulation trajectory. This tool calculates the three components of the binding energies: (1) molecular mechanics, (2) polar solvation, and (3) nonpolar. The complete trajectory of 100 ns was split into 5 ns intervals that resulted in 20 frames. These 20 frames were used in binding energy calculation. Figure S3 shows the MMPBSA binding energy during the simulation on 20 frames. The total binding energy for the complex (protein-9-6) is always negative, thus confirming the stability of the complex. The molecular mechanics component contributed most toward producing negative binding energy that negates the polar solvation energy, as shown in Figure S3 (Table ).

Table 2

MMPBSA Energy Components of 3CL Protease Binding with both the Ligands (9 and 6)a

	protease <−> (’9′ + ’6′)		9′ <−>’6’
energy components	energy value (kJ/mole)	standard deviation (+/−)	energy value (kJ/mole)	standard deviation (+/−)
van der Waals energy	–245.802	7.25	–0.014	0.001
electrostatic energy	24.404	0.974	0.023	0.017
polar solvation energy	150.551	5.215	1.043	1.333
SASA (solvent accessible surface area) energy	–26.525	0.722	–0.039	0.032
binding energy	–97.355	7.167	1.04	1.316

The binding energy between both the ligands is also shown in the table.

Average values for each energy component and the total binding energy are shown in Table S1 between the protein and the ligands (9 and 6). Additionally, the average binding energies between these two ligands are also shown in Table S1. The standard deviation shown indicates that binding energy does not highly fluctuate during the simulation time; the average MMPBSA binding energy for the protein-9-6 complex is −97.35 kJ/mol with ±7.1 standard deviation. The MMPBSA binding energy between 9 and 6 is not negative due to the ∼15 Å distance but has a very small magnitude of 1.31 kJ/mol, indicating that both the ligands are thermodynamically stable in their respective binding site.

Binding energy contribution from each residue of the protein in binding both the ligand was also calculated to determine the importance of residue. Figure S4 shows the average energy contribution of each residue in binding. Active site residues are marked in orange to show the energetic behavior of these residues in binding the ligands. GLN189 from the active site has the minimum most energy, −5.6 kJ, which showed its strong interaction with ligand 9. Only one active site residue, ARG188, showed positive energy in the binding. Catalytic dyads are marked red in Figure S4. Both HIS41 and CYS145 showed negative average binding energy across the 100 ns simulation. CYS145 stands at second rank with −5.2 kJ average binding energy and thus contributes significantly to protein–ligand interaction.

Dose–Response Assay

For up to seven concentrations ranging from 0.1 to 100 μM, we found a dose-dependent inhibition of 3Cl protease activity (Figure a). Compounds 9 and 6 resulted in a half-maximum effective concentration (EC50) of 45 μM. Within the measured spectrum, there was no observable cytotoxicity, with the concentration causing 50% cytotoxicity (CC50) being more than 100 μM (Figure ).

Figure 5

(a) Efficacy of the compound was calculated by the percentage inhibition of SARS-CoV-2 3CL protease by the specified concentrations of compounds 6 and 9 in combination ranged from 0.1 to 100 μM. Data are shown as mean ± S.D. for three technical replicates. Binding kinetics of (b) compound 6, (c) compound 9, and (d) combination of compounds 6 and 9 to the 3CL Protease. Kinetic analysis of 3CL Protease binding was performed by injecting different known concentrations of aggregates (from 0.31 to 5 μM) onto the immobilized carboxymethyl dextran-coated CM5 sensor chip was used with amine-coupling chemistry. All measurements were performed at 25 °C with a flow rate of 30 μL/min using HBS-EP buffer with association time 60 s, followed by 60 s dissociation phase. Kinetic constants (KDs) were calculated from the sensorgrams using the 1:1 fit model using BIA Evaluation 2.0.1 (Cytiva) software. The blue line indicates the actual curve, and the black line indicates the fitted curve of the sensorgram.

Figure 6

Cytotoxicity study of the compounds 9 and 6 individually and in combination at 1, 10, 50, and 100 μM concentrations.

SPR Assay

SPR assay indicates the binding of two molecules. Here, compounds 9 and 6 were allowed to bind to the protein individually and then simultaneously, as we showed using enzyme-based assay. The smaller the KD value, the greater the binding affinity of the ligand for its target. The binding kinetics of both the compounds were determined at different concentrations. The association and dissociation curves of the SPR assay are shown in Figure b–d for compound 6, compound 9, and compound 9 and 6, respectively. It is observed that both compounds 9 and 6 exhibit the same order of the kinetic affinities individually (high KD value), thereby indicating a relatively lower affinity to the target, as shown in Table . However, the KD value decreases significantly when the two compounds are tested together, confirming their synergistic effect (Table ).

Table 3

In Vitro Target Binding of Compounds 9 and 6 and their Mixture to the 3CL Proteasea

sample	ka	kd	KD (M)	Chi2 (RU2)	U-value
9	1.39 × 105	0.01643	1.18 × 10–7	0.186	1
6	1.96 × 104	0.0123	6.32 × 10–7	2.24	1
9 and 6	621,100	0.04057	6.53 × 10–8	0.192	3

Estimated Kd values were accessed with Biacore X100TM by Cytiva.

Cell Viability

The MTT assay was used to evaluate the cytotoxicity of the compounds used in this study. In a 96-well plate, 10,000 cells per well were seeded overnight at 37 °C and 5% CO2. The next day, the wells were incubated in triplicate with the required amount of compound for a period of 36 h at 37 °C and 5% CO2. The media was substituted with a 100 μL MTT solution of 0.5 mg/mL prepared in complete media and incubated at 37 °C for 2 h. The supernatant was discarded, and 150 μL of dimethyl sulfoxide (DMSO) was used to dissolve the formed crystals in each well for 15 min. The absorbance was calculated at 595 nm in an ELISA plate reader. Cell viability was recorded for the individual and combined application of compounds 9 and 6, as shown in Figure . At 1, 10, 50, and 100 μM concentrations, compound 9 has cell viabilities of more than 94, 93, 97, and 92%, respectively, as shown in Figure . Compound 6 showed cell viability that was greater than 93, 84, 77, and 68% at 1, 10, 50, and 100 μM concentrations, respectively. However, when compounds 9 and 6 were used in combination, cell viability was more than 91, 82, 73, and 62% at 1, 10, 50, and 100 μM concentrations, respectively.

Methods

General

Compounds used in this study are obtained from a commercial supplier; characterization reports of these compounds are included in Table S1. Their purification sheet is provided from the vendor; all compounds have a reported purity >98% as determined by high-performance liquid chromatography.

Collection of the Protein Structure

The three-dimensional structure of the protein was collected from the protein data bank.[38] The 3CL protease structure was solved and deposited in the PDB with ID: 6M2N20. The PDB structure of protease (6M2N) deposited in the database has 2.2 Å resolution as per X-ray crystallography. The experimental structure co-ordinates used in this study are crystallized complexes with a novel inhibitor molecule targeting their binding site. Shuanghuanglian, a Chinese patent medicine, has two critical bioactive components, baicalein and baicalin. This structure has baicalein docked with 3CL protease. Protease has three domains, d-I (residues 8–101), d-II (residues 102–184), and d-III (residues 201–303).[39] Domains I and II are the substrate-binding domains; baicalein was productively ensconced in the core of the substrate-binding pocket by interacting with two catalytic residues, the crucial S1/S2 subsites.

Sourcing of Small Molecules

Molecules that are used against 3CL protease in this study were collected from three different resources. These are primarily the repository for natural origin compounds (1) ZINC natural molecules (https://zinc15.docking.org/substances/subsets/natural-products/), (2) BIOFACQUIM (https://biofacquim.herokuapp.com/), and (3) Chemfaces (http://www.chemfaces.com/). Molecules from BIOFACQUIM and Chemfaces were docked directly, while ZINC molecules were first screened using pharmacophores of the bound baicalein structure, and then top hits were docked to identify potential hits.

Virtual Screening

Ligand-based pharmacophores were generated using the Pharmagist server[40] using various inhibitors or binders suggested for COVID-19. The top 50 compounds from the pharmacophore screening were docked using AutoDock Vina to determine their binding energies. However, compounds from BIOFACQUIM and Chemfaces were directly docked using AutoDock Vina because of the small population, and top hits were selected based on the binding energy.

Molecular Docking

Molecular docking between small compounds and 3CL protease was performed using AutoDock Vina and Autodock-4.[41,42] Protein and ligand were prepared using AutoDock Tool and Marvin sketch; multiple docking was performed using Raccon application of AutoDock. Hydrogen atoms were added to the system, and Gasteiger charges were assigned to prepare the ligand. The pKa value of the ligand was determined using the ChemAxon tool integrated with Marvin sketch application. Docking used the genetic algorithm to explore the binding space. The number of individuals generated in each cycle of the genetic algorithm was 150, while a maximum of 2,500,000 energy evaluations were performed. Docking program used a maximum of 27,000 generations, where only the top individual survives to the next generation with a 0.02 rate of mutation and a 0.8 rate of crossover.

MD simulation of the selected docked complex was performed using GROMACS 4.6.2[43] having CHARMM27 force field.[44] Simulations were carried out under physiological conditions. Topology and parameters for small molecules compatible with the CHARMM all-atom force field were generated using the CGenFF program.[45] Hydrogen atoms were added to the RdRp protein under physiological conditions and placed in a solvated box at a 1.4 nm distance from the wall. The protein–ligand solvated complex was energetically minimized using 5000 steps of the steepest deepest method. A time step of 2 fs was used using the SHAKE algorithm. Constant temperature (NVT) and pressure ensemble (NPT) conditions were applied to the system for 100 ps and 1 ns, respectively, to attain the equilibrium state. Ligand and protein molecules were constrained during the equilibrium phase. After the equilibrium phase, a 100 ns all-atom simulation was performed using V-rescale temperature coupling[46] for an external heat bath with a 0.1 ps time constant for the protein and ligand, while pressure coupling was performed using Parrinello–Rahman with a time constant of 2 ps; long-range electrostatic was dealt with PME (particle mesh ewald) method.

Enzyme Assay

A SARS-CoV-2 3CL protease assay kit was purchased from BPS Biosciences (San Diego, CA, USA), which includes recombinant 3CL protease enzyme, 3CL protease substrate, 3CL protease assay buffer, GC376 as an inhibitor control, and a black, low binding microtiter plate for assay performance. Briefly, 3CL protease was thawed on ice and diluted to a concentration of 5 ng/μL in assay buffer containing 1 mM DTT. Thirty microliters of diluted 3CL protease enzyme solution (150 ng) was added to the microtiter well plate in triplicate for positive, inhibitor, and test samples. Inhibitor compounds were dissolved in DMSO at a 100-fold higher concentration than the actual concentration test in the assay. Ten microliters of the test inhibitor was added in triplicate in the required amount to the designated well. GC376 at 100 mM was added to the inhibitor control in triplicate. Ten microliters of 5% DMSO was added to the blank and positive control wells, and the plate was incubated for 30 min at room temperature with slow shaking. The substrate for 3CL protease was diluted in assay buffer at a concentration of 250 μM. The reaction was started by adding 10 μL of the substrate solution to each well. The final concentration of the 3CL protease substrate in a 50 μL reaction was 50 μM. The plate was sealed with a plate sealer and incubated at room temperature overnight. The fluorescence intensity was measured using SpectraMax microplate readers with Ex = 360 nm/Em = 460 nm. The blank value was subtracted from all other values when analyzing the final data.

The binding kinetics of 3CL Protease with the targeted molecules and their combination were evaluated by using a Biacore X-100 system with CM5 chips (Cytiva). The 3CL Protease protein was immobilized on the chip with a concentration of 33 μg/mL (diluted by 0.1 mM NaAc, pH 4.5). For all measurements, the same running buffer was used, which consisted of 20 mM HEPES, pH 7.5, 150 mM NaCl, and 0.005% tween-20 with association time 60 s, followed by a 60s dissociation phase. A blank channel of the chip was used as negative control. Serially diluted protein solutions (0.31, 0.62, 1.25, 2.5, and 5 μM) then flowed through the chip surface. The multicycle binding kinetics was analyzed with the Biacore X-100 Evaluation Software and fitted with a 1:1 binding model.

Conclusions

In this paper, we present a case of mutual synergistic effect of two natural origin chemical compounds that act simultaneously on 3CL protease of SARS-CoV-2. Initially, the individual efficacy of both compounds was evaluated, and both of them displayed low inhibitory activity against the 3CL protease. Docking and simulation suggested that binding of the molecules is short-term, resulting in their lower efficacy. Furthermore, MD simulation explained that the removal of ligands from the binding site of the protein is prevented when both the compounds are docked with the protein molecule. These compounds are docked at different binding sites. Compound 9 has flavonoid scaffold and is bound at the active site of the 3CL protease protein and exhibits strong interactions with the catalytic residues (HIS41 and CYS145). However, compound 6 has the naphthalene chemical scaffold and sits 15 Å away from the catalytic site. The binding of 6 prevents the removal of 9 from the catalytic site, making it a longer-term interaction. The combined assay was performed initially at 100 and 50 μM concentrations of each compound. At both these concentrations, inhibition activity percentage was high at 97 and 85%. Later, the multiconcentration experiment of the combined effect was performed from 0.1 to 100 μM of each compound, and it was observed that even a 0.5 μM concentration causes ∼20 to 25% inhibition activity. SPR assay was also performed to find the binding affinity of these compounds both individually and in combination. KD values obtained in SPR suggest 10-fold stronger binding of the combination.

Overall, our study showcases successful application of molecular simulation for identifying potential inhibitors for the 3CL protease activity of SARS-CoV-2. The synergetic effect of two compounds predicted by the in silico study is successfully validated by results from the enzyme and SPR assay.