Design and Evaluation of Andrographolide Analogues as SARS-CoV-2 Main Protease Inhibitors: Molecular Modeling and in vitro Studies

Introduction

The coronavirus disease 2019 (COVID-19), caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has had a profound global impact, resulting in over seven million deaths worldwide (https://covid19.who.int/, accessed on March 11, 2025). While immunity can be achieved through natural infection or vaccination—both of which help reduce mortality risks—the emergence of new variants presents significant challenges. These variants often exhibit increased transmissibility, infectivity, and the potential to evade existing immunity.^1,2 Consequently, there is an urgent need to discover and develop effective antiviral agents to slow disease progression and complement ongoing vaccine development efforts.

SARS-CoV-2 is a positive-sense, single-stranded RNA virus belonging to the Coronaviridae family and the Betacoronavirus genus. Its genome encodes nine accessory proteins and four structural proteins: the spike (S), envelope (E), membrane (M), and nucleocapsid (N) proteins. In addition, the RNA genome of SARS-CoV-2 includes sixteen non-structural proteins (NSP1 to NSP16) expressed as part of the replicase polyproteins 1a and 1ab (pp1a and pp1ab).³ Key non-structural proteins—such as the main protease (M^pro, also known as 3-chymotrypsin-like protease or 3CL^pro), papain-like protease (PL^pro), helicase (NSP13), and RNA-dependent RNA polymerase (RdRp)—are potential targets for COVID-19 therapies.⁴ Specifically, M^pro plays a crucial role in viral replication by cleaving pp1a and pp1ab at multiple sites to produce essential mature NSP required for forming the RNA replicase-transcriptase complex.⁵ Furthermore, the absence of closely related homologs in humans enhances M^pro’s attractiveness as a target for COVID-19 treatment, minimizing the risk of off-target effects.⁶

Structurally, SARS-CoV-2 M^pro is composed of three distinct domains: domain I (residues 8–101), domain II (residues 102–184), and domain III (residues 201–303), with domains II and III linked by a long loop.⁶ These domains form a homodimer, a configuration that exhibits significantly higher hydrolytic activity than the monomeric form.⁷ A key factor in the formation of this dimeric structure is the salt bridge interaction between Arg4 of one monomer and Glu290 of the other.⁸ Additionally, the catalytic dyad, consisting of His41 and Cys145 located in the respective domains I and II, acts as a nucleophile and a general acid/base, playing a crucial role in the cleavage of the viral polyprotein.⁷

In response to the COVID-19 pandemic, substantial efforts have been made to develop drugs specifically targeting SARS-CoV-2 M^pro. Among these, nirmatrelvir, a potent covalent inhibitor of M^pro, combined with ritonavir under the brand name Paxlovid, has been approved in several countries for COVID-19 treatment.⁹ Clinical studies show that nirmatrelvir can reduce hospitalization risk by 89% in unvaccinated COVID-19 patients.¹⁰ Another effective M^pro inhibitor, ensitrelvir (S-217622, marketed as Xocova), is a non-covalent, non-peptidomimetic drug that has shown clinical benefits in reducing COVID-19 severity and has been approved for use in Japan.¹¹ In addition, leritrelvir¹² and simnotrelvir¹³ have received conditional approval in China. While the success of nirmatrelvir and ensitrelvir demonstrates the effectiveness of targeting M^pro in combating SARS-CoV-2, concerns remain about the potential emergence of drug resistance due to mutations in the viral genome. Mutations such as S144A, E166V, and A173V in SARS-CoV-2 M^pro variants have been shown to reduce the potency of both nirmatrelvir and ensitrelvir.¹⁴ Therefore, developing new antivirals that can effectively inhibit SARS-CoV-2 M^pro is a pressing need.

Natural products offer a promising candidate for drug discovery due to their structural diversity and biological activity.¹⁵ The semisynthetic modification of natural products allows for the introduction of various functional groups onto their core structures, enhancing their pharmacological properties.¹⁶ Andrographolide, the primary active compound in Andrographis paniculata, is known for its broad-spectrum antiviral activities, including against SARS-CoV-2.¹⁷ In our previous in vitro study on semisynthetic derivatives of andrographolide, we observed promising outcomes. Compounds 3k, 3l, 3m, and 3t demonstrated more than 50% inhibitory activity against SARS-CoV-2 M^pro at a concentration of 10 μM.¹⁸ Furthermore, molecular dynamics (MD) simulations combined with binding free energy calculations showed that the (S)-enantiomer of compound 3l exhibited the strongest binding affinity to M^pro. This suggests potential for further structural optimization to enhance binding efficiency and bioactivity. Specifically, the γ-crotonolactone, benzyl, and p-fluorobenzyl groups in (S)-3l were well-positioned within the S1, S2, and S4 pockets of M^pro, respectively, which may contribute to its potent inhibitory activity.

Building on these promising results, we employed a structure-based drug design approach to develop new andrographolide analogues using compound (S)-3l as the template. Drawing inspiration from ensitrelvir—the first oral non-covalent SARS-CoV-2 M^pro inhibitor featuring a 2,4,5-trifluorobenzene pharmacophore—we designed novel analogues by incorporating a 2,4,5-trifluorobenzene group at the P2 position (Figure 1). This moiety is known to occupy the hydrophobic S2 pocket of SARS-CoV-2 M^pro, mimicking the P2 position of ensitrelvir and potentially enhancing binding affinity.¹⁹ To further explore the effects on protein–ligand interactions, we also varied the chemical substituents at the P2 position. Therefore, the aim of this study is to design and evaluate C-12 dithiocarbamate andrographolide derivatives as potential inhibitors of SARS-CoV-2 M^pro through a combined in silico and in vitro approach. We hypothesize that the introduction of a dithiocarbamate and a 2,4,5-trifluorobenzene moiety will enhance binding affinity and inhibitory activity against SARS-CoV-2 M^pro by stabilizing key interactions within the active site. The findings provide a basis for developing potent andrographolide-based analogues, contributing to ongoing antiviral drug discovery efforts against COVID-19.

Figure 1 Rational design of novel andrographolide analogues based on ensitrelvir framework. The atomic labels of the core structure of andrographolide are also provided.

Materials and Methods

Computational Details

Preparation of Protein and Ligand

The previously docked structure of the andrographolide analogue (S)-3l in complex with SARS-CoV-2 Mpro ¹⁸ served as the template for our structure-based drug design approach. To create models of the newly designed andrographolide analogues (1–6) bound to protomer A (chain A) of SARS-CoV-2 M^pro, we performed structural modifications on the p-fluorobenzene moiety of the original (S)-3l molecule. These modifications transformed (S)-3l into the corresponding analogues 1–6 (Figure 2), utilizing the small molecules tool available in Discovery Studio Visualizer (BIOVIA, San Diego, CA, USA). For comparison, the structure of ensitrelvir bound to SARS-CoV-2 M^pro was retrieved from the Protein Data Bank (PDB ID: 7VU6). The protonation states of ionizable amino acids (Asp, Glu, Lys, Arg, and His) were considered at pH 7.4 and set based on our previous study.^18,20 The chemical structures of 1–6 were constructed using the Gaussview 5.0 program and were structurally optimized at the B3LYP/6-31G(d) level of theory using the Gaussian 16 program.²¹ Next, the electrostatic potential (ESP) charges of each ligand were calculated through single-point calculations at the HF/6-31G(d) level of theory. The antechamber and parmchk2 modules implemented in AMBER22 were used to generate the restrained ESP (RESP) charges and the missing molecular parameters for the ligands, respectively. The AMBER ff19SB force field ²² and the generalized AMBER force field 2 (GAFF2)²³ were adopted to treat the bonded and nonbonded parameters for the protein and ligand, respectively. Missing hydrogen atoms were added to the protein–ligand complex using the LEaP module of AMBER22. Each system was subsequently solvated with the OPC explicit solvation model,²⁴ ensuring a minimum buffer thickness of 10 Å, and was electrostatically neutralized with sodium (Na⁺) counterions. The newly added hydrogen atoms in the protein–ligand complex and solvent molecules underwent energy minimization through 1000 steps of steepest descent (SD), followed by 2000 steps of conjugate gradient (CG) methods. Next, the protein and ligand were energy-minimized using the same approach, with 1000 steps of SD and 2000 steps of CG, while the solvent molecules were kept fixed. Finally, the entire system underwent complete energy minimization, employing 1000 iterations for the SD method and 2000 iterations for the CG method.

Figure 2 Chemical structures of designed andrographolide analogues 1–6.

Molecular Dynamics Simulations

To comprehensively assess the stability and dynamic behavior of the protein–ligand complexes in an aqueous environment, MD simulations were conducted using the SANDER and PMEMD modules of the AMBER22 software package.²⁵ The modeled system was executed under periodic boundary conditions with the isothermal–isobaric (NPT) scheme, as detailed in our previous study.^18,26 Long-range electrostatic interactions were addressed using the particle-mesh Ewald summation method,²⁷ while a 10-Å cutoff distance was established for nonbonded interactions. The SHAKE algorithm²⁸ was employed to maintain all hydrogen atom bonds at a constant length, allowing for a 2-fs time step of integration. To control the pressure and temperature of each system, the Berendsen barostat,²⁹ with a pressure-relaxation time of 1 ps, and the Langevin thermostat,³⁰ with a damping frequency of 2 ps⁻¹, were applied. The simulated system was gradually thermalized from 10 to 300 K over 200 ps using the canonical ensemble (NVT), with positional restraints of 30.0 kcal/mol·Å² on the C_α atoms of the protein. Following this, the complex underwent NPT equilibration for 1300 ps, which included four steps of restrained MD simulations. During these steps, the restraints on the atoms of the protein active site were progressively decreased from 30 to 20, 10, and finally 5 kcal/mol·Å², followed by an additional 500 ps without any restraints. Subsequently, the entire system was subjected to NPT ensemble MD simulations at 300 K and 1 atm until a total simulation time of 500 ns was reached. MD trajectories were collected every 10 ps for analysis. The post-dynamic trajectories were examined for their structural features using the CPPTRAJ module³¹ of AMBER22. Per-residue decomposition free energy calculations were performed using the MMPBSA.py module³² to further quantify the contributions of individual residues to ligand binding, allowing for a more detailed understanding of the key interactions driving inhibitory activity. 3D visualizations were generated using the UCSF Chimera³³ and ChimeraX³⁴ programs.

Experimental Details

General Information for the Synthesis of Compound 1

All chemicals were purchased from commercial suppliers and used without further purification. Proton nuclear magnetic resonance (¹H NMR) spectra were recorded on a Bruker Avance 400 MHz spectrometer using chloroform-d (CDCl₃) as the solvent. Chemical shifts (δ) are reported in parts per million (ppm) relative to tetramethylsilane (TMS, δ 0.00) or CDCl₃ (δ 7.26) as internal standards. Carbon nuclear magnetic resonance (¹³C NMR) spectra were recorded on a Bruker Avance 100 MHz spectrometer using CDCl₃ as the solvent. Chemical shifts (δ) are reported in ppm relative to CDCl₃ (δ 77.0) as the internal standard. High-resolution mass spectrometry (HRMS) data were obtained at Naresuan University, Thailand. Analytical thin-layer chromatography (TLC) was performed on precoated silica gel 60F-254 plates (E. Merck, Darmstadt, Germany). Column chromatography used silica gel 60 PF254 (E. Merck, Darmstadt, Germany). Melting points were determined using a Griffin melting point apparatus and are reported uncorrected.

Preparation of 14-Acetyl-3, 19-Isopropylidene Andrographolide as a Precursor for the Synthesis of Compound 1

The synthetic pathway for producing 14-acetyl-3,19-isopropylidene andrographolide, which serves as a precursor for the synthesis of compound 1 is shown in Scheme 1.

Scheme 1 The procedure for the synthesis of 14-acetyl-3,19-isopropylidene andrographolide.

Step 1: Protection of Andrographolide with 2,2-Dimethoxypropane

Andrographolide (2.0016 g, 5.71 mmol) was dissolved in acetone (60 mL). 2,2-Dimethoxypropane (3.5 mL) and PPTS (0.7177 g) were added, and the reaction mixture was stirred at room temperature for 2 hours. The progress of the reaction was monitored by TLC using 70% ethyl acetate (EtOAc)/hexane as the eluent. Upon completion, the reaction was quenched by adding saturated sodium hydrogen carbonate (NaHCO₃) solution until the pH of the aqueous layer was neutral. The mixture was extracted with EtOAc (x3) and the combined organic layers were washed with saturated sodium chloride (NaCl) solution, dried over anhydrous sodium sulfate (Na₂SO₄), and concentrated under reduced pressure to afford the crude product. The crude product was purified by column chromatography using 50% EtOAc/hexane as the eluent.³⁵ The 3,19-isopropylidene andrographolide was obtained as a white solid in 91% yields (2.03 g) (R_f = 0.78) and was further performed acetylation reaction.

3,19-isopropylideneandrographolide; ¹H NMR (400 MHz, CDCl₃): δ6.96 (t, J = 7.0 hz, 1H), 5.03 (brt, 1H), 4.91 (s, 1H), 4.62 (s, 1H), 4.45 (dd, J = 10.0, 6.0 hz, 1H), 4.27 (dd, J = 10.0, 1.5 hz, 1H), 3.96 (d, J = 11.5 hz, 1H), 3.49 (dd, J = 9.0, 4.0 hz, 1H), 3.18 (d, J = 11.5 hz, 1H), 2.60–2.54 (m, 2H), 2.42 (dt, J = 12.5, 3.0 hz, 2H), 2.05–1.92 (m, 1H), 1.85–1.78 (m, 1H), 1.78–1.70 (m, 2H), 1.41 (s, 3H), 1.36 (s, 3H), 1.34–1.22 (m, 4H), 1.20 (s, 3H), 0.96 (s, 3H).

Step 2: Acetylation

The 3,19-isopropylidene andrographolide (0.5007 g, 1.28 mmol) was dissolved in EtOAc (3.5 mL) in a round-bottom flask. Triethylamine (0.8047 mL, 5.77 mmol) in EtOAc (4.0 mL) was added, followed by dropwise addition of acetyl chloride (320 μL, 4.49 mmol). The reaction mixture was stirred at room temperature for 30 minutes. The progress of the reaction was monitored by TLC using 50% EtOAc/hexane as the eluent. Upon completion, the reaction was quenched by the addition of saturated NaHCO₃ solution until the aqueous layer reached neutral pH. The mixture was extracted with EtOAc (x3) and the combined organic layers were washed with saturated NaCl solution, dried over Na₂SO₄, and concentrated under reduced pressure to afford the crude product. The crude product was purified by column chromatography using 20% EtOAc/hexane as the eluent. Each collected fraction was analyzed by TLC. The desired fractions were combined and concentrated under reduced pressure, then further dried under vacuum to obtain product as yellow liquid in 89% yields (0.49 g) (R_f=0.92). 14-acetyl-3,19-isopropylidene andrographolide; ¹H NMR (400 MHz, CDCl₃): δ 7.03 (td, J = 6.9, 1.7 hz, 1H), 5.93 (d, J = 6.1 hz, 1H), 4.89 (s, 1H), 4.55 (dd, J = 11.3, 6.1 hz, 1H), 4.54 (s, 1H), 4.24 (dd, J = 11.2, 2.0 hz, 1H), 3.95 (d, J = 11.6 hz, 1H), 3.50 (dd, J = 8.6, 4.0 hz, 1H), 3.18 (d, J = 11.6 hz, 1H), 2.47–2.37 (m, 3H), 2.12 (s, 3H), 2.05–1.93 (m, 2H), 1.88–1.68 (m, 3H), 1.40 (s, 3H), 1.36 (s, 3H), 1.32–1.22 (m, 4H), 1.19 (s, 3H), 0.94 (s, 3H).

Synthesis and Characterization of Andrographolide Analogue (Compound 1)

The synthetic pathway for the production of compound 1 is outlined in Scheme 2 as follows:

Scheme 2 The procedure for the synthesis of compound 1.

Step 1: Dithiocarbamate Generation

A solution of 1-(2,4,5-trifluorobenzyl)piperazine (370 mg, 1.61 mmol) in CH₃CN (3 mL) was prepared at room temperature. Carbon disulfide (CS₂) was then added at 3 °C. The mixture was stirred at room temperature until complete conversion was confirmed by TLC.³⁶

Reaction with Andrographolide Derivative

To the reaction mixture, 14-acetyl-3,19-isopropylidene andrographolide (695 mg, 1.61 mmol) was added and stirred at room temperature. After confirming complete conversion by TLC, the reaction mixture was cooled, and excess silica gel was added, followed by stirring overnight.

Isolation of Product

The reaction mixture was filtered, washed with EtOAc, and the organic solvent was removed under reduced pressure to yield the compound 12-((1-(2,4,5-trifluorobenzyl)piperazine-4-carbonothioyl)thio)-3,19-isopropylidene-14-deoxyandrographolide as a pale brown solid in 80% yield (873 mg, 1.28 mmol). This compound was then used in the second step to synthesize compound 1.

Precursor of Compound 1

12-((1-(2,4,5-trifluorobenzyl)piperazine-4-carbonothioyl)thio)-3,19-isopropylidene-14-deoxyandrographolide: Pale brown solid (873 mg, 80%); ¹H NMR (400 MHz, CDCl₃): δ 7.42 (t, J = 1.4 hz, 1 h), 7.25–7.19 (m, 1H), 6.95–6.88 (m, 1H), 5.21 (dd, J = 12.3, 4.2 hz, 1H), 4.98 (s, 1H), 4.80 (dd, J = 4.1, 1.5 hz, 2H), 4.75 (s, 1H), 4.30 (brs, 2H), 3.94 (d, J = 11.6 hz, 1H), 3.92 (brs, 2H), 3.54 (s, 2H), 3.44 (dd, J = 9.2, 4.0 hz, 1H), 3.15 (d, J = 11.6 hz, 1H), 2.53 (s, 4H), 2.41 (dt, J = 12.7, 2.4 hz, 1H), 1.98–1.93 (m, 2H), 1.92–1.68 (m, 6H), 1.40 (s, 3H), 1.35 (s, 3H), 1.31–1.21 (m, 3H), 1.18 (s, 3H), 0.89 (s, 3H); ¹³C NMR (100 MHz, CDCl₃): δ 195.83, 172.63, 157.76, 155.32, 148.62, 147.72, 146.65, 134.63, 132.02, 119.49, 108.96, 99.29, 70.41, 64.14, 54.51, 53.98, 52.97, 52.21, 46.40, 39.02, 38.70, 38.40, 38.12, 34.54, 27.79, 25.69, 25.49, 23.79, 16.53; HRMS (ESI) m/z calcd for C₃₅H₄₅F₃N₂O₄S₂ [M+H]⁺ 679.2851, found 679.2854.

Step 2

A solution of 12-((1-(2,4,5-trifluorobenzyl)piperazine-4-carbonothioyl)thio)-3,19-isopropylidene-14-deoxyandrographolide (347 mg, 0.51 mmol) in 2 M HCl (1 mL) was stirred for 15 minutes. After complete conversion was confirmed by TLC, the mixture was filtered, washed with EtOAc, and the organic solvent was removed under reduced pressure. Compound 1 (12-((1-(2,4,5-trifluorobenzyl)piperazine-4-carbonothioyl)thio)-3,19-hydroxy-14-deoxyan drographolide) was obtained as a white solid in 67% yield (219 mg, 0.34 mmol).

Compound 1

(12-((1-(2,4,5-trifluorobenzyl)piperazine-4-carbonothioyl)thio)-3,19-hydroxy-14-deo xyandrographolide: White solid (219 mg, 67%); ¹H NMR (400 MHz, CDCl₃): δ 7.54–7.39 (m, 2H), 7.06 (s, 1H), 5.10 (d, J = 9.6 hz, 1H), 4.96 (s, 1H), 4.84 (s, 2H), 4.69 (s, 1H), 4.32–4.17 (m, 5H), 3.44–3.29 (m, 4H), 2.56–2.06 (m, 7H), 1.74–1.40 (m, 6H), 1.25 (s, 3H), 1.05–0.82 (m, 3H), 0.64 (s, 3H); ¹³C NMR (100 MHz, CDCl₃): δ 196.96, 172.44, 158.19, 148.70, 148.44, 145.95, 131.27, 122.53, 122.35, 108.41, 106.41, 80.30, 70.27, 64.13, 64.04, 55.27, 54.10, 53.32, 46.87, 39.10, 38.85, 36.59, 26.66, 26.46, 23.93, 22.27, 15.50; HRMS (ESI) m/z calcd for C₃₂H₄₁F₃N₂O₄S₂ [M+H]⁺ 639.2538, found 639.2539.

SARS-CoV-2 M^pro in vitro Assay

The assays for protease activity and inhibition were conducted following the methods outlined in previous studies.^20,37,38 In brief, enzyme kinetics were carried out using SARS-CoV-2 M^pro at a concentration of 0.2 µM and the fluorogenic peptide substrate E(EDANS)TSAVLQSGFRK(DABCYL) (Biomatik), with excitation and emission wavelengths set at 340 and 490 nm, respectively. During the initial screening for inhibitory activity, enzymatic activity was assessed both with and without compound concentrations of 100 μM. To determine the half-maximal inhibitory concentration (IC₅₀), the initial rate of substrate cleavage (25 μM) was measured at various concentrations of the selected compounds. The IC₅₀ values were analyzed using GraphPad Prism 9 (GraphPad Software, Inc., San Diego, CA, USA). The Cheng-Prusoff equation³⁹ was applied to calculate the inhibitory constant (K_i) using a previously reported K_m value of 51 μM.³⁷

Cytotoxicity

HepG2 (ATCC HB-8065), HEK293T (ATCC CRL-3216), and HaCaT (Catalog No. T0020001) cell lines, obtained from the American Type Culture Collection (ATCC, USA) and AddexBio (San Diego, USA), respectively, were seeded in 96-well clear plates (Corning, Acton, MA, USA) at densities of 1×10⁴, 1.5×10⁴, and 1×10⁴ cells per well. The cells were cultured for 24 hours at 37°C with 5% CO₂ in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibco, Detroit, MI, USA) supplemented with 10% fetal bovine serum (FBS, Merck, Schuchardt, Darmstadt, Germany) and 1% penicillin-streptomycin (P/S, Gibco). For HEK293T cells, 2 mM L-glutamine (Hyclone, GE Healthcare, USA) was also added. After incubation, cells were treated with compound 1 at concentrations of 100, 50, 25, 12.5, 6.25, 3.13, and 1.56 μM for 24 hours under the same conditions. Ensitrelvir (Sigma-Aldrich, St. Louis, MO, USA) served as the positive control. Following treatment, the medium was replaced with serum-free medium containing MTT reagent [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide] (Abcam, Cambridge, UK) at a final concentration of 0.5 mg/mL. The cells were incubated for an additional 3 hours at 37°C with 5% CO₂. The MTT-containing medium was then discarded, and DMSO (Merck, Schuchardt, Darmstadt, Germany) was added to dissolve the formazan crystals. Absorbance was measured at 570 nm using a microplate reader (Biometrics Technologies, Wilmington, DE, USA). CC₅₀ values were calculated using GraphPad Prism 9 (GraphPad Software, Inc., San Diego, CA, USA).

Results and Discussion

This study employed a combined in silico and in vitro approach to design and evaluate C-12 dithiocarbamate andrographolide analogues as potential SARS-CoV-2 M^pro inhibitors. Computational analysis guided the selection of structural modifications that enhance ligand binding affinity and stability within the M^pro active site. Among the tested compounds, compound 1 emerged as the most promising candidate, demonstrating notable inhibitory activity against SARS-CoV-2 M^pro while maintaining a favorable cytotoxicity profile. These findings highlight the potential of andrographolide derivatives for antiviral development and underscore the need for further studies to optimize their potency and safety for clinical applications, as discussed in the following sections.

Molecular Dynamics of Designed Compounds Bound to SARS-CoV-2 M^pro

Structural Stability and Compactness of Protein–Ligand Complex

The stability of the protein–ligand complexes was assessed using root-mean-square deviation (RMSD), which is commonly employed to evaluate the dynamic stability of modeled complexes in an aqueous environment.^40–44 As shown in Figure 3A, all simulated systems exhibited similar RMSD trends, characterized by increased fluctuations during the initial 100 ns, followed by an equilibration phase between 400 and 500 ns. Note that all modeled complexes displayed slightly lower RMSD values during the simulation time compared to ensitrelvir, suggesting that these complexes could form stable and equilibrate more rapidly in the aqueous environment.

Figure 3 Analysis of structural dynamics. (A) Time evolution of the backbone root-mean-square deviation (RMSD) of the amino acids within 5 Å of the ligand. (B) The calculated radius of gyration (Rg) of the designed compounds and ensitrelvir in complex with SARS-CoV-2 Mpro.

The compactness of the protein structure in complex with the designed compounds was further assessed by calculating the radius of gyration (R_g) of the protein C_α atoms, as depicted in Figure 3B. The average R_g values ranged from 25.33 to 25.45 Å, indicating minimal variation across all modeled complexes. This consistency indicates that the protein maintained a tightly packed structure throughout the simulations. Additionally, the R_g values are in agreement with previous MD studies on andrographolide analogues¹⁸ and nirmatrelvir,⁴⁵ further supporting the stability of these complexes. Based on these findings, the structural coordinates from the final 50 ns (450–500 ns) of MD simulations were selected for subsequent analyses.

Atomic Contacts and Intermolecular Hydrogen Bonds of Protein–Ligand Complex

To briefly investigate the binding capability of the designed compounds toward SARS-CoV-2 M^pro, the number of atomic contacts (# contacts) between each compound and the active site residues was calculated using a 3.5 Å cutoff. Any atoms within this distance were counted as contacts. As shown in Figure 4, the # contacts for these complexes remained relatively stable throughout the simulation period. During the final 50 ns, the # contacts for compounds 1–6 ranged from 3 to 20, with ensitrelvir exhibiting the highest # contacts at 27. Notably, compound 1 exhibited similar # contacts (20 ± 5 atoms) to ensitrelvir (27 ± 6 atoms), suggesting that compound 1 may bind more favorably to the active site residues of SARS-CoV-2 M^pro compared to the other designed compounds. Additionally, the # contacts of compound 1 are comparable to other reported M^pro inhibitors, such as lapatinib⁴⁶ and fragmented-lapatinib aminoquinazoline analogues,⁴⁷ further supporting its potential as a promising inhibitor.

Figure 4 Number of atomic contacts (# contacts) between the ligand and the active site residues plotted over the entire 500-ns simulation period. The average # contacts from the final 50 ns is represented as mean ± SD.

Next, we performed an analysis of hydrogen bond (H-bond) formation to quantitatively evaluate the intermolecular interactions between the active site amino acids and the designed compounds. The analysis was based on two geometric criteria: (i) an acceptor-donor distance ≤ 3.5 Å and (ii) an acceptor-H-donor angle ≥ 120°. As illustrated in Figure 5, the designed andrographolide analogues demonstrated limited H-bonding with the active site residues of SARS-CoV-2 M^pro, forming only 1 to 3 H-bonds. In comparison, ensitrelvir established 5 to 6 H-bonds. These findings could provide valuable insights into the mechanisms of compound recognition and the energetic contributions to ligand binding, which are discussed further in the following section.

Figure 5 Time evolution of the number of intermolecular hydrogen bonds (H-bonds) between the surrounding amino acids and the ligand.

Water Accessibility at the SARS-CoV-2 M^pro Active Site

To assess water accessibility at the SARS-CoV-2 M^pro active site, solvent-accessible surface areas (SASAs) were calculated for amino acid residues within 5 Å of each ligand. Only protomer A contained the ligand, while protomer B represented the apo form (Figure 6A). The time evolution of SASA for each system is shown in Figure 6B. The results show that the average SASAs for the apo form of compounds 1–6 ranged from approximately 890 to 1200 Ų, while ensitrelvir showed an average SASA of 1086 ± 101 Ų. Upon ligand binding in the active site, the average SASAs for all simulated models decreased by approximately 110 to 414 Ų. This observation is consistent with previous MD studies on lopinavir/ritonavir,²⁰ peptidomimetic inhibitors (including compounds N3, 11a, 13b, and 14b),⁴⁸ and halogenated baicalein,⁴⁹ all of which reported significant reductions in SASA during the binding process. Among the designed compounds, 1, 2, and 5 showed slightly higher SASA values than ensitrelvir but remained lower than the other analogues. These findings indicate that the designed compounds are well-embedded within the SARS-CoV-2 M^pro binding pocket.

Figure 6 (A) Representation of the SARS-CoV-2 Mpro homodimer, showing protomer A (light blue) with the ligand bound and protomer B (light pink) in the unbound state. (B) Calculated solvent-accessible surface area (SASA) of the six modeled complexes during 500-ns MD simulation. Amino acids within a 5-Å radius of the ligand were selected for SASA calculations. The average SASA values (mean ± SD) for protomer A and protomer B, derived from the last 50 ns of individual MD simulations, are depicted in light blue and pink, respectively.

End-State Binding Free Energy Calculations

To evaluate the binding strength of the designed compounds against SARS-CoV-2 M^pro in a dynamic system, 500 snapshots from the last 50 ns of MD simulations were analyzed to calculate ∆G_bind using the solvated interaction energy (SIE) method. The SIE method is an end-point, physics-based scoring function that directly calibrates solvation terms for binding energy calculations in protein–ligand complexes.^50,51 Compared to the commonly used MM/PB(GB)SA method, SIE provides a faster and more efficient estimation of binding free energies by eliminating the need for computationally intensive entropy calculations.

∆G_bind is a widely used thermodynamic parameter for evaluating protein–ligand binding affinity, making it a valuable criterion for identifying promising drug candidates.^52,53 As shown in Table 1, the ΔG_bind values for the designed compounds ranged from −9.69 to −6.73 kcal/mol, with ensitrelvir exhibiting a ΔG_bind of −9.92 kcal/mol. Notably, the predicted ΔG_bind for ensitrelvir closely aligned with its experimentally derived value of −10.37 kcal/mol, calculated from its IC₅₀ (0.049 ± 0.001 µM⁵⁴) using the formula . This agreement could validate the predictive method’s reliability and support the accuracy of the results. In addition, van der Waals interactions (ΔE_vdW) and molecular surface-based energies (ΔG_cavity) played a key role in ligand binding, significantly contributing to the overall binding affinity. Although coulombic interactions (ΔE_coul) were favorable, their contribution was often offset by unfavorable reaction field energies (ΔG_RF), reducing their overall impact on binding energy. For further investigation, compound 1, which exhibited the lowest ΔG_bind, along with the highest # contacts (see Figure 4) among the designed compounds, was selected for experimental evaluation of its biological activity against SARS-CoV-2 M^pro, as detailed in the following section. Besides, the calculated energy terms in Table 1 highlight the critical role of specific noncovalent interactions in drug recognition. The designed compounds showed significantly greater contributions from van der Waals interactions compared to electrostatic interactions (8- to 22-fold lower), consistent with previous studies emphasizing the hydrophobic nature of the SARS-CoV-2 M^pro binding pocket.^20,48,49 In addition, we can draw conclusions about the impact of P2 substituents on binding affinity against SARS-CoV-2 M^pro. Among the designed compounds (1–6), those containing isobutyl (2) and cyclohexyl (3) moieties exhibited weaker binding affinity to M^pro, suggesting that aliphatic and bulkier aliphatic groups may not favorably interact with the S2 pocket of M^pro. In contrast, the introduction of benzyl and fluorobenzyl moieties enhanced binding affinity, following the trend: 2,4,5-trifluorobenzyl (1) > benzyl (4) > p-fluorobenzyl (6) > m-fluorobenzyl (5). These findings highlight the critical role of benzyl and fluorobenzyl substitutions in enhancing the binding affinity of andrographolide analogues toward SARS-CoV-2 M^pro.

Table 1 ΔGbind Values and Their Energetic Components (in kcal/mol) for the Complexes of the Designed Compounds and Ensitrelvir with SARS-CoV-2 Mpro Estimated Using the End-State SIE Method

It should be noted that a few reports have investigated andrographolide derivatives as potential COVID-19 therapeutics, exploring different mechanisms of action. While our study focuses on direct M^pro inhibition, other research highlights alternative antiviral pathways. For instance, Schulte et al 2022,⁵⁵ demonstrated that andrographolide derivatives can modulate the KEAP1/NRF2 signaling pathway, enhancing cellular defense mechanisms and thereby reducing viral replication. Meanwhile, Thomas et al 2024,⁵⁶ employed a cheminformatics-driven approach to identify andrographolide derivatives as dual inhibitors of SARS-CoV-2 methyltransferases (nsp14 and nsp16). These enzymes play a key role in viral RNA capping, a process essential for evading host immune recognition.

In vitro Evaluation of Anti-SARS-CoV-2 M^pro Activity and Predicted Key Binding Residues

The inhibitory activity of compound 1 against SARS-CoV-2 M^pro was evaluated using an M^pro inhibition assay with a fluorogenic substrate. Initially, the relative inhibition at a concentration of 100 μM was measured for both andrographolide and compound 1. Compound 1 demonstrated approximately 70% inhibition of M^pro activity (Figure 7A), which is more potent than the parent andrographolide and previously reported inhibitors, such as caffeic acid ester derivatives (which showed 58.0–68.8% relative inhibition at 100 μM).⁴⁸ This finding suggests that the designed compound exhibits enhanced binding affinity and inhibitory activity against SARS-CoV-2 M^pro. Subsequently, a dose-response curve was generated to determine the IC₅₀ value of compound 1. As depicted in Figure 7B, the IC₅₀ value of compound 1 was calculated to be 72.0 ± 3.9 μM (K_i ~ 48.3 μM), placing it within a similar micromolar range as known potent non-covalent inhibitors, including X77 and ML188,⁵⁷ as well as semi-synthetic compounds such as baicalein and brominated baicalein (IC₅₀ values of 63 ± 1 and 56 ± 1 μM, respectively).⁴⁹ These results highlight the potential of compound 1 as a promising scaffold for further development as a SARS-CoV-2 M^pro inhibitor.

Figure 7 In vitro evaluation of compound 1’s inhibitory activity and analysis of key amino acid residues in SARS-CoV-2 Mpro involved in its binding. (A) Percentage of relative inhibition of SARS-CoV-2 Mpro activity in the presence of 100 μM compound 1 compared to its parent compound, andrographolide. (B) Dose-response curve of Relative Fluorescence Units (RFU) generated during the enzyme’s initial rate period over time (%RFU/s), used to calculate the IC50 value of compound 1 against SARS-CoV-2 Mpro. (C) Per-residue free energy decomposition of key binding residues, presented with a surface representation indicating their contribution levels (blue to green).

To better understand the mechanism of inhibitory action, the amino acids acting as hot-spot residues were predicted using the MM/GBSA method. Residues with a free energy decomposition value of ≤ −1.0 kcal/mol were identified as key binding amino acids. As illustrated in Figure 7C, six amino acids—H41, T45, M49, M165, D187, and R188—were considered critical contributors to the binding of compound 1. Among these, M49 and M165 were the most influential, consistent with previously reported key residues involved in the binding of 12-dithiocarbamate-14-deoxyandrographolide analogues,¹⁸ brominated baicalein,⁴⁹ and a newly designed compound MI-23.⁵⁸ These findings indicate that compound 1 establishes strong intermolecular interactions with the catalytic residue H41, as well as with the substrate-binding residues, which are crucial for M^pro inhibition.⁵⁹

Furthermore, the inhibitory mechanism of compound 1 can be attributed to its ability to form several non-covalent interactions within the active site of SARS-CoV-2 M^pro. The 2,4,5-trifluorobenzene moiety in compound 1 enhances hydrophobic interactions, particularly with M49 and M165, which are crucial for ligand stabilization. Furthermore, H-bonding interactions involving H41 further strengthen binding affinity and enzyme inhibition. These key interactions align with those observed in other non-covalent M^pro inhibitors, such as ensitrelvir,¹⁹ α-mangostin derivatives,⁶⁰ and ebsulfur derivatives.⁶¹ Given that compound 1 exhibits comparable binding affinity to known inhibitors, it is likely that its mechanism of action involves interference with the substrate-binding region, thereby preventing the proteolytic processing of viral polyproteins.

A structural comparison between compound 1 and ensitrelvir reveals both similarities and differences in their binding modes and inhibitory activities. Ensitrelvir is a non-covalent SARS-CoV-2 M^pro inhibitor that incorporates a 2,4,5-trifluorobenzyl moiety, the same as in compound 1, which enhances hydrophobic interactions within the S2 pocket of M^pro. Both compounds also interact with Met49, a key residue involved in ligand stabilization. However, ensitrelvir exhibits significantly greater potency, with an IC₅₀ of 0.013 μM, compared to compound 1 (IC₅₀ = 72.0 μM). This difference is likely due to the optimized P1 and P1’ substitutions of ensitrelvir, which include a 1-methyl-1H-1,2,4-triazole moiety in the S1 pocket, enhancing H-bonding with H163, and a 6-chloro-2-methyl-2H-indazole moiety in the P1’ position, which contributes to additional hydrophobic interactions.¹⁹ While compound 1 exhibits promising inhibition of SARS-CoV-2 M^pro, further structural optimization, such as modifying the P1 moiety to strengthen interactions with H163 and optimizing a P1’ moiety that forms stronger hydrophobic contacts with surrounding residues could help bridge the potency gap between compound 1 and clinically approved inhibitors.

Cytotoxic Effects on Representative Cell Lines

Understanding the cytotoxicity profile of the identified hit compound is important for assessing its therapeutic potential during drug development. In this study, HaCaT cells (human keratinocyte cells), HEK293T cells (human embryonic kidney cells), and HepG2 cells (human liver cancer cells) were used to represent three major organ systems. As shown in Table 2 and Figures S1–S3 in the supplementary data, compound 1 demonstrated moderate to low cytotoxicity, with CC₅₀ values of 13.24 ± 2.25 µM in HaCaT cells, 41.02 ± 2.62 µM in HEK293T cells, and 42.26 ± 1.16 µM in HepG2 cells. While compound 1 exhibited higher toxicity than ensitrelvir, its cytotoxicity levels were comparable to those of its parent compound, andrographolide, across all tested cell lines. These results indicate that compound 1 enhances binding to SARS-CoV-2 M^pro without significantly compromising the viability of liver and kidney cells, which are essential for drug metabolism and excretion. However, the increased sensitivity of HaCaT cells to compound 1 (CC₅₀ = 13.24 ± 2.25 µM) may limit its potential use in dermatological applications. Compared to other M^pro inhibitors, compound 1 exhibits a toxicity profile similar to several non-covalent M^pro inhibitors, including α-mangostin analogues,⁶⁰ sulfonamide chalcones,⁶² and flavonoid derivatives.⁶³ While its cytotoxicity profile appears favorable, further selectivity studies against human proteases and long-term toxicity evaluations in preclinical models are essential to confirm its safety for therapeutic development.

Table 2 Cytotoxic Effects of Compound 1 Indicated by CC50 Values (µM) in Three Different Cell Lines: HaCaT Cells (Human Keratinocyte Cells), HEK293T Cells (Human Embryonic Kidney Cells), and HepG2 Cells (Human Liver Cancer Cells)

However, this study has certain limitations that should be mentioned. While the in vitro enzyme inhibition assay provides valuable insights into the inhibitory activity of compound 1 against SARS-CoV-2 M^pro, it does not fully capture its antiviral efficacy in a cellular environment. Another significant limitation is the lack of cytotoxicity testing in lung cell lines, which are more physiologically relevant for SARS-CoV-2 infection. Although cytotoxicity assessments in HaCaT, HEK293T, and HepG2 cells provide initial safety insights, they may not accurately reflect the effects of the compound in lung tissues, the primary site of SARS-CoV-2 infection. Future studies should incorporate lung-derived epithelial cells to provide a more relevant and comprehensive assessment of cytotoxicity. Moreover, the absence of in vivo validation limits the ability to assess the pharmacokinetic and toxicity profiles, which are critical for therapeutic development. To address these limitations, further investigations, including cell-based antiviral assays and animal model studies, are necessary to confirm its clinical applicability and therapeutic potential.

Conclusions and Future Perspectives

This study highlights the successful design and evaluation of C-12 dithiocarbamate andrographolide analogues as SARS-CoV-2 M^pro inhibitors, supporting the development of novel antiviral therapeutics. Using a structure-based drug design approach, six analogues (1–6) were developed and analyzed through MD simulations, which revealed their stable binding within the M^pro active site. Among them, compound 1, incorporating a 2,4,5-trifluorobenzene moiety, exhibited the most favorable binding energy (ΔG_bind = −9.69 kcal/mol) and formed strong interactions with key active site residues (H41, M49, and M165). Experimental validation further supported the computational findings, as compound 1 demonstrated ~70% inhibition of M^pro at 100 μM with an IC₅₀ value of 72 μM. Cytotoxicity assessments indicated moderate to low toxicity, comparable to andrographolide, suggesting a favorable safety profile in vitro. These findings underscore the potential of andrographolide derivatives as promising scaffolds for SARS-CoV-2 antiviral drug development. To enhance its therapeutic potential, future studies should focus on optimizing compound 1 by improving its potency, selectivity, and pharmacokinetic properties through structural modifications and analogue development. In addition, further in vitro and in vivo evaluations are necessary, such as testing against clinically relevant M^pro mutants, conducting antiviral assays in SARS-CoV-2-infected cells, and assessing pharmacokinetic and toxicity profiles in animal models. These investigations will provide critical insights into its therapeutic feasibility and potential for clinical translation.

Beyond its relevance to SARS-CoV-2, the findings of this study have broader implications for antiviral drug discovery. The structure-based optimization of andrographolide derivatives could serve as a basis for developing inhibitors against other coronaviruses, such as SARS-CoV-1 and Middle East Respiratory Syndrome Coronavirus (MERS-CoV), which share structural and functional similarities in their M^pro.^64,65 Additionally, the computational drug design strategies utilized in this study may also be applicable to other viruses that depend on cysteine proteases for their life cycle, such as human rhinovirus, enterovirus 71 and picornavirus.^66–68 Given the ongoing risk of emerging coronaviruses, the identification of natural product-derived inhibitors as broad-spectrum antiviral agents presents a promising strategy for pandemic preparedness. Future research could extend these findings by evaluating optimized compound 1 and its analogues against a range of viral proteases to assess their potential for broader antiviral applications. This study contributes to the ongoing efforts to develop effective antiviral therapeutics against COVID-19 and emerging coronaviruses, reinforcing the significance of natural product-derived inhibitors in antiviral drug discovery campaigns.