Multi-ligand molecular docking, simulation, free energy calculations and wavelet analysis of the synergistic effects between natural compounds baicalein and cubebin for the inhibition of the main protease of SARS-CoV-2

Abstract

Combination drugs have been used for several diseases for many years since they produce better therapeutic effects. However, it is still a challenge to discover candidates to form a combination drug. This study aimed to investigate whether using a comprehensive in silico approach to identify novel combination drugs from a Chinese herbal formula is an appropriate and creative strategy. We, therefore, used Toujie Quwen Granules for the main protease (M^pro) of SARS-CoV-2 as an example. We first used molecular docking to identify molecular components of the formula which may inhibit M^pro. Baicalein (HQA004) is the most favorable inhibitory ligand. We also identified a ligand from the other component, cubebin (CHA008), which may act to support the proposed HQA004 inhibitor. Molecular dynamics simulations were then performed to further elucidate the possible mechanism of inhibition by HQA004 and synergistic bioactivity conferred by CHA008. HQA004 bound strongly at the active site and that CHA008 enhanced the contacts between HQA004 and M^pro. However, CHA008 also dynamically interacted at multiple sites, and continued to enhance the stability of HQA004 despite diffusion to a distant site. We proposed that HQA004 acted as a possible inhibitor, and CHA008 served to enhance its effects via allosteric effects at two sites. Additionally, our novel wavelet analysis showed that as a result of CHA008 binding, the dynamics and structure of M^pro were observed to have more subtle changes, demonstrating that the inter-residue contacts within M^pro were disrupted by the synergistic ligand. This work highlighted the molecular mechanism of synergistic effects between different herbs as a result of allosteric crosstalk between two ligands at a protein target, as well as revealed that using the multi-ligand molecular docking, simulation, free energy calculations and wavelet analysis to discover novel combination drugs from a Chinese herbal remedy is an innovative pathway.

Keywords: ADME/T, absorption, distribution, metabolism, excretion and toxicity; COVID-19; COVID-19, Coronavirus disease 2019; Combination drug therapy; Computer simulation; Computers molecular; H-bonds, hydrogen bonds; LD50, median lethal dose; MD, molecular dynamics; MM-PBSA, molecular mechanics Poisson Boltzmann surface area; Mpro, main protease; Natural products; PAINS, Pan-assay interference compounds; RCO, inter-residue contact order; RMSF, root-mean-square-fluctuation; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; SMILES, Simplified Molecular-input Line-entry System; TCMSP, traditional Chinese medicine systems pharmacology database and analysis platform; TQG, Toujie Quwen Granule; Virus diseases.

Graphical abstract

Fig. 1

Crystal structures of SARS-CoV-2 main protease in complex with compounds from the herbal remedy Toujie Quwen Granules. a) Docking pose interactions between all compounds and the protease. The yellow wire mesh is the surface area of the active binding sites of the protease, whereas the surface of primary binding pockets is labeled in blue wire mesh. His41 and Cys145 are the key active binding site residues. b-h) interactions formed between the seven potential inhibitors identified from single ligand molecular docking and surrounding residues. The yellow wire mesh is the surface area of their binding pockets. Hydrogen bonds are presented by green dashed lines and their distances are labeled respectively. Corresponding compound names were presented in Supplementary Table S3.

Fig. 2

Synergistic binding modes of the four compounds directly connected to HQA004 (baicalein) against the main protease of SARS-CoV-2. HQA004 is presented in green color sticks. a) Overview of the structure of HQA004 bound to the four compounds. b-e) 3D/2D binding poses of the four compounds complex with HQA004 and their surrounding residues that were predicted by multiple ligand molecular docking. Hydrogen bonds are presented by green dashed lines while hydrophobic bonds are in pink. Corresponding compound names were presented in Supplementary Table S3.

Fig. 3

Simulated diffusion pathway of synergistic ligand CHA008. a) Multi-frame overlays for compound CHA008 from the initial synergistic docking position (red licorice representation) to the final stably-bound position (blue licorice representation), shown every 2 ns. M^pro is shown as light grey ribbons, with compound HQA004 as yellow spheres. b-d) Time series plots for minimum-contact distances between CHA008 and Ser1, Pro108, and Asn203 respectively. Graphical insets below the plots illustrate ligand positions. ‘Checkpoint’ residues on M^pro are shown as multicolored spheres. HQA004 is shown in yellow, and CHA008 is shown in light green. e-f) 2D ligand interaction diagrams are displayed for CHA008 bound at the initial docking-predicted site (e) and the final potential allosteric site (f).

Fig. 4

CHA008 promotes the conformational stability of HQA004. a-b) Multi-frame overlays for compound HQA004 from the initial docking pose (red licorice representation) to the final MD-simulated pose (blue licorice representation), shown every 2 ns. M^pro is shown as light grey ribbons, with compound CHA008 as light green spheres. Singly-bound HQA004-M^pro complex is shown in a), and doubly-bound HQA004-CHA008-M^pro complex is shown in b). c-d) Structurally-aligned multi-frame overlays for the HQA004 molecule within the c) singly-bound and d) doubly-bound M^pro complex. e-f) Time-series plots of the central C–C dihedral angle of HQA004 for e) singly-bound and f) doubly-bound M^pro complex.

Fig. 5

CHA008 promotes the binding affinity of HQA004. a) Time-series plot of the total number of inter-atomic contacts between HQA004 and M^pro for singly- (black) and doubly-bound M^pro with CHA008 (grey) or LQA003 (red). b) Per-residue binding free energy contribution of M^pro residues towards HQA004, calculated using MM-PBSA, for singly-bound (black bars) and doubly-bound M^pro with CHA008 (grey bars) or LQA003 (red bars). c-d) Graphical illustrations of the impact of CHA008 on specific contacts between HQA004 and the surrounding active site residues for the c) singly-bound and d) doubly-bound M^pro complex. HQA004 is shown in multicolored licorice format, while active site residues are shown as spheres-and-sticks, and color-coded according to residue name. 2D ligand interaction diagrams are inset underneath.

Fig. 6

Proposed allosteric mechanism of CHA008. a) Per-residue differences in root-mean-square-fluctuation (ΔRMSF) for doubly-bound M^pro relative to singly-bound M^pro, with negative values indicating enhanced stability due to CHA008 synergistic binding, and vice versa. b) M^pro shown in ribbon representation, with segments color-coded according to ΔRMSF value using a scheme comprising red (lowest ΔRMSF), white (intermediate ΔRMSF), and blue (highest ΔRMSF). HQA004 is shown as yellow spheres, with CHA008 as green spheres. c) Multi-frame overlays of the Cα trace of the first principal component (vibrational mode) of the doubly-bound M^pro complex, color-coded using a red-white-blue scheme to represent concerted motions.

Fig. 7

Wavelet power spectra for the ligand-bound chain A of Mpro, calculated for △RMSF (a) and △RCO (b), where the differences were calculated at each residue by subtracting the RMSF/RCO value obtained for the singly-bound (CHA-free) from that of the doubly-bound (CHA bound) complex. The black contours show a 10% significance level, and the color code reflects the strength of the spectrum, ranging from blue (low power) to red (high power). Regions of special interest described in the manuscript text are indicated by the dashed box.

Fig. 8

Synergistic mechanism network model between HQA004 and other 206 compounds from Toujie Quwen Granules. For corresponding herbal names, please refer to Table 1; Corresponding compound names were presented in Supplementary Table S2.