Discovery of C-12 dithiocarbamate andrographolide analogues as inhibitors of SARS-CoV-2 main protease: In vitro and in silico studies.

A global crisis of coronavirus disease 2019 (COVID-19) pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has impacted millions of people's lives throughout the world. In parallel to vaccine development, identifying potential antiviral agents against SARS-CoV-2 has become an urgent need to combat COVID-19. One of the most attractive drug targets for discovering anti-SARS-CoV-2 agents is the main protease (Mpro), which plays a pivotal role in the viral life cycle. This study aimed to elucidate a series of twenty-one 12-dithiocarbamate-14-deoxyandrographolide analogues as SARS-CoV-2 Mpro inhibitors using in vitro and in silico studies. These compounds were initially screened for the inhibitory activity toward SARS-CoV-2 Mpro by in vitro enzyme-based assay. We found that compounds 3 k, 3 l, 3 m and 3 t showed promising inhibitory activity against SARS-CoV-2 Mpro with >50% inhibition at 10 μM. Afterward, the binding mode of each compound in the active site of SARS-CoV-2 Mpro was explored by molecular docking. The optimum docked complexes were then chosen and subjected to molecular dynamic (MD) simulations. The MD results suggested that all studied complexes were stable along the simulation time, and most of the compounds could fit well with the SARS-CoV-2 Mpro active site, particularly at S1, S2 and S4 subsites. The per-residue decomposition free energy calculations indicated that the hot-spot residues essential for ligand binding were T25, H41, C44, S46, M49, C145, H163, M165, E166, L167, D187, R188, Q189 and T190. Therefore, the obtained information from the combined experimental and computational techniques could lead to further optimization of more specific and potent andrographolide analogues toward SARS-CoV-2 Mpro.

Introduction

The coronavirus disease 2019 (COVID-19) pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread across the globe with >6 million deaths and confirmed infection cases of over 470 million as of March 22, 2022 (https://covid19.who.int/). Since the outbreak of SARS-CoV in 2002 and MERS-CoV (Middle East respiratory syndrome) in 2012 [1], the emergence of the highly pathogenic SARS-CoV-2 has totally changed the scenario worldwide, which has never been seen in the past century. This causative agent of COVID-19 mainly attacks the respiratory system and potentially the other parts of the body (e.g., gastrointestinal tract, renal system and nervous system) [2], [3], [4]. Some infected patients could develop severe pneumonia and acute respiratory distress syndrome with trouble breathing, which necessarily requires ventilator assistance and intensive care, especially in the elderly, immunocompromised, and those with underlying conditions [5]. To date, many vaccines have been developed by both academics and pharmaceutical companies and then authorized for emergency use [6]. Although there have been continuous efforts to develop COVID-19 vaccines, there are still some issues that could affect their practical uses such as vaccine inequity in low- and middle-income countries and vaccine hesitancy. More importantly, there is concern with respect to the effectiveness of vaccines in preventing viral infection with emerging new variants of concern [7]. After the alpha, beta, gamma and delta SARS-CoV-2 variants, the highly contagious omicron (B.1.1.529) variant with a large number of mutations and some deletions was emerged and reported to World Health Organization (WHO) from South Africa, and this variant became spread very quickly to many countries worldwide [8]. Thus, in parallel to COVID-19 vaccine development, antiviral drug discovery for the treatment of SARS-CoV-2 infection is still urgently needed.

SARS-CoV-2 is an enveloped, positive-sense, single-stranded RNA betacoronavirus of the family Coronaviridae as SARS-CoV and MERS-CoV. A large genome of SARS-CoV-2 consists of approximately 30 kb in length, where its entire genome is highly similar to SARS-CoV with a sequence identity of ∼80% [9]. One‐third of the genome encodes for several accessory proteins and four structural proteins comprising a spike glycoprotein (S), an envelope protein (E), a membrane protein (M) and a nucleoprotein (N), which play a crucial role in virion assembly [10]. Whereas two overlapping open reading frames (ORFs), namely ORF1a and ORF1b, encode the respective replicase polyproteins, so-called pp1a and pp1ab, that locate at the two‐thirds of genomic RNA [11]. Both pp1a and pp1ab are further proteolytically cleaved into 16 mature non-structural proteins (NSPs) by two viral proteases, papain‐like cysteine protease (PLpro) and main protease (Mpro), encoded in ORF1a. Consequently, these liberated NSPs assemble to the viral RNA polymerase complex, which is required for genome transcription and replication of new virion. In particular, SARS-CoV-2 Mpro or known as 3C‐like protease (3CLpro) or NSP5 with the molecular weight of 33.8 kDa is the pivotal enzyme of virus, which is responsible for cleaving polyproteins at least 11 out of the 14 conserved scissile junctions. Owing to the indispensable function in the viral life cycle and none of the human homologous proteins, Mpro has been gained much attention as a promising therapeutic target for COVID-19 treatment [12], [13]. Currently, several approved drugs, drug candidates, pharmacologically active compounds [14] and fragment-based design [15] have been widely explored by structure-based virtual and high-throughput screening aimed at slowing down or stopping SARS-CoV-2 growth in the host cells via Mpro inhibition.

In general, a homodimeric SARS-CoV-2 Mpro composes of three domains (domains I-III) linking domains II and III through a long loop [16], as shown in Fig. 1A. The native substrate‐binding site contains the H41 and C145 catalytic dyad located between domains I and II. Note that dimerization is required to form the catalytically active Mpro and trigger protease activity [17]. The Mpro homodimer is stabilized mainly through a salt bridge interaction between two domain III and maintained by two NH2 termini (N-fingers: residues 1–7) compressed between two protomers [14]. To date, numerous co-crystallized ligands and compound fragments have been deposited in the Protein Data Bank (PDB) and clarified to block the catalytic site of SARS-CoV-2 Mpro [14], [15], [18], [19].

Fig. 1

(A) 3D structure of SARS-CoV-2 Mpro in complex with non-covalent inhibitor X77 (PDB ID: 6W63) at the enzyme binding pocket of protomer A (green). Domains I, II and III of protomer B are represented by red, cyan and blue, respectively. (B) 2D chemical structure of andrographolide with its atomic labels. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Among all promising candidate compounds, phytochemicals have played an essential role in new drug discovery and development for decades [20]. The ent-labdane andrographolide (Fig. 1B) isolated from the “King of Bitters” plant Andrographis paniculata is one of the prominent natural products that has been attracted attention in recent years. It has been shown to exert a wide range of pharmacological activities, including anti-inflammatory, antioxidant, anticancer, hepatoprotective and particularly antiviral effects [21], [22], [23], [24]. A. paniculata extract is traditionally used to alleviate common cold symptoms, fever and diarrhea, as well as upper respiratory and inflammatory diseases [25]. The clinical trials show that there are no serious adverse effects observed in patients [26], [27]. More interestingly, andrographolide was computationally predicted to bind to the specific targets regarding SARS-CoV-2 infection, including the host angiotensin-converting enzyme 2 (ACE2) and three viral proteins (S protein, Mpro and RNA-dependent RNA polymerase) [28], [29], [30], [31]. Recently, in vitro enzyme-based assay identified that andrographolide exhibited an inhibitory effect toward SARS-CoV-2 Mpro [32]. Further study reported that A. paniculata extract and pure andrographolide showed significant anti-SARS-CoV-2 activity against infected human lung epithelial cell (Calu-3). Besides, the favorable cytotoxicity profiles among major organ cell models indicated the high-​potential development of both A. paniculata extract and andrographolide to fight SARS-CoV-2 infection either as a monotherapy or in combination regimens [33]. However, andrographolide indeed possesses multi-targeting property (i.e., one single molecule binding to multiple targets) [34], which might limit the specific application for COVID-19 treatment. To enhance and optimize the potency and selectivity of andrographolide, chemical modification of the parent compound has been typically conducted [35]. Specifically, a series of 12-dithiocarbamate-14-deoxyandrographolide analogues [36] (Fig. S1 and Table S1, Supplementary Materials) is of particular interest due to the extension of the long chain at C-12 of andrographolide (Fig. 1B). We hypothesized that this introduction would make the stronger binding interactions to the S2, S3 and S4 subsites of SARS-CoV-2 Mpro, as previously guided by in silico molecular docking [29]. It is noteworthy that this series of andrographolide analogues was successfully synthesized in one pot under mild conditions and formerly tested anticancer activity against several cancer cell lines [36]. Herein, to extend the anti-SARS-CoV-2 Mpro potential of such analogues, we firstly screened the promising compounds using in vitro SARS-CoV-2 Mpro inhibition assay. Then, the binding mode and interaction of the selected compounds toward SARS-CoV-2 Mpro were studied by computational methods. Ultimately, we hope that these findings can be helpful for the future design and development of andrographolide analogues in the fight against COVID-19.

Materials and methods

Experimental details

SARS-CoV-2 Mproin vitro assay

As previously described, the enzyme-based assay for SARS-CoV-2 Mpro was carried out [37]. SARS-CoV-2 Mpro was expressed and purified using a method as previously reported for SARS-CoV Mpro [38]. The enzyme kinetic was conducted with SARS-CoV-2 Mpro at 0.2 µM using the fluorogenic substrate E(EDANS)TSAVLQSGFRK(DABCYL) (Biomatik), in which the excitation and emission wavelength were measured at 340 and 490 nm, respectively. For the initial screening of inhibitory activity, the enzymatic activity was measured in the presence and absence of compound concentration at 10 and 100 μM. The initial rate in the absence of an inhibitor was used for normalization. For the half-maximal inhibitory concentration (IC50) determination, the initial rate of substrate (25 μM) cleavage was measured when the selected compounds were present at various concentrations. The IC50 value was fitted with GraphPad Prism 8 (GraphPad Software, Inc., La Jolla, CA, USA). The corresponding inhibitory constant (Ki) value was calculated using the Cheng-Prusoff equation [39] with the previously reported Km value (51 μM) [37].

Computational details

Preparation of protein and ligand

The starting X-ray structure with 2.1 Å resolution of homodimeric SARS-CoV-2 Mpro bound to the co-crystallized non-covalent inhibitor (X77) was obtained using PDB ID 6W63 [40]. Since all the 12-dithiocarbamate-14-deoxyandrographolide analogues were expected to be non-covalent inhibitors, thus, the 6W63 structure was chosen because there were a limited number of crystal structures of SARS-CoV-2 Mpro in complex with non-covalent inhibitors at the time we performed the computational studies. The 6W63 structure also showed relatively high resolution and clear electron density of ligand and the active site residues among the available crystal structures. Besides, X77 is a non-covalent inhibitor that occupies all four main pockets (S1′, S1, S2 and S4) of the SARS-CoV-2 Mpro active site, thereby allowing other ligands to bind to the substrate-binding cleft of the enzyme. The protonation states of all ionizable residues (K, R, D, E and H) were assigned in accordance with our previous study [41]. All basic residues (K and R) were set as the positively charged form, whereas all acidic residues (D and E) were set as the negatively charged form. The H residue was assigned based on the prediction from the H++ web server [42] at pH 7.4. This was except for the catalytic residue H41, which was set as the neutral form with protonated δ-NH (HID type of AMBER format) according to the mechanistic consideration of cysteine proteases [43]. The chemical structures of 12-dithiocarbamate andrographolide analogue candidates (Fig. 2) were built utilizing Gaussview 5.0 program. These ligands were firstly optimized at the B3LYP level of theory with 6-31G* basis set using Gaussian09 program [44]. Afterward, the partial atomic charges and empirical force field parameters for each optimized ligand were developed according to our standard procedure [45], [46] as follows. The electrostatic potential (ESP) charges for each of the optimized geometries were computed by single-point calculation at the HF/6-31G* level of theory. Then, the restrained ESP (RESP) charges were derived by converting the ESP charges using the antechamber module of AMBER20 [47]. The parmchk2 module was utilized to generate the missing molecular parameters of the ligand based on the general AMBER force field version 2 (GAFF2) [48]. Whereas the AMBER ff19SB force field [49] was adopted to treat the bonded and non-bonded parameters of the protein.

Fig. 2

Chemical structures of focused 12-dithiocarbamate-14-deoxyandrographolide analogues, where the labels of heteroatoms are also indicated. Note that a chiral configuration at the chlorcyclizine scaffold of 3l and 1-phenyltetrahydroisoquinoline moiety of 3t was presented and referred to as (R,S)-3l and (R,S)-3t for further discussion.

Molecular docking

To predict the mode of binding between SARS-CoV-2 Mpro and each andrographolide analogue, molecular docking using Autodock Vina [50] with rigid docking was carried out. Firstly, the optimized structures of andrographolide analogues (see detail above) were viewed and converted to MOL2 format using Gaussview 5.0 program. The MOL2 files were then changed to pdbqt format using AutoDock Tools version 1.5.6 [51]. For the preparation of a receptor model in docking simulations, all solvent molecules and the inhibitor in the 6W63 protein were removed. The structure of X77 (CID: 145998279) was downloaded from the PubChem web server and subjected to optimization at the B3LYP/6-31G* level using the Gaussian09 program, which was further used for the redocking experiment. The cubical grid box of 30×30×30 Å at a spacing of 1 Å was defined using the co-crystallized inhibitor X77 as the centroid by keeping other defaults setup. The exhaustiveness value for docking was set to 64. Note that the redocking simulation was performed to validate the docking protocol. The superimposition between the docked pose and the co-crystallized ligand showed that the docked X77 with the highest binding affinity (i.e., the first pose with the most negative value) was well aligned with the crystal conformation (root-mean-square deviation, RMSD of 0.99 Å, see Fig. S2). Thus, the docking method with these settings could reproduce the crystal binding pose of the ligand in the protein active site and was dependable to apply for such a system. Further, each focused compound was separately docked into the active site of SARS-CoV-2 Mpro using the same docking procedure. The optimum pose between SARS-CoV-2 Mpro and each ligand was chosen for subsequent molecular dynamics (MD) simulation. Apart from considering the docking score, it should be mentioned that the optimum docking model for each candidate was selected using two more specific criteria as follows: (i) showing γ-crotonolactone ring in the Mpro S1 pocket and (ii) displaying aryl ring in the Mpro S2 pocket. Otherwise, if the docking simulations fail to generate any poses aligning with these criteria, ligand will be manually positioned within the Mpro active site to maximize the interactions between the preferential scaffold of compound and the Mpro S1 and S2 pockets (discussed later in section 3.2 of results and discussion). This is because the S1 and S2 pockets of SARS-CoV-2 Mpro have been found to be steadily occupied by γ-lactam moiety (structurally similar to γ-crotonolactone ring) and the bulkier group such as cyclohexyl and fluorobenzyl moiety of the reported peptidomimetic inhibitors, respectively [18], [19].

Molecular dynamics simulations

The addition of missing hydrogen atoms of the protein–ligand complex was performed using the LEaP module implemented in AMBER20 according to the aforementioned molecular parameters. Each simulated system was solvated in a simulation box of the OPC explicit solvation model [52] with a minimum buffer thickness of 10 Å, containing ∼19,000 water molecules. The sodium counterions (8 Na+) were randomly added to keep the whole system neutral. To remove the bad contacts, the added hydrogen atoms and water molecules were minimized using 1000 iterations of steepest descent (SD) minimization continued by 2000 iterations of conjugated gradient (CG) with the remaining atoms fixed. Afterward, the protein and ligand were subjected to energy minimization using the same approach of SD (1000 steps) and CG (2000 steps) with the solvent molecules fixed. In the final step, the whole complex was fully energy-minimized with 1000 steps and 2000 steps of SD and CG methods, respectively.

MD simulations of the protein–ligand complexes were simulated with the modules SANDER and PMEMD of the AMBER20 software package [53]. The modeled system was run under the periodic boundary condition with the isothermal–isobaric (NPT) scheme, as described in our previous study [41]. In brief, the long-range electrostatic interactions were treated using the particle mesh Ewald’s summation method [54], whereas a 10-Å cutoff distance was set for nonbonded interactions. All the hydrogen atoms bonds were maintained at a constant bond length with the SHAKE algorithm [55], allowing 2 fs time step of integration. The Berendsen barostat [56] with a pressure-relaxation time of 1 ps and the Langevin thermostat [57] with a damping frequency of 2 ps−1 were applied to control the pressure and temperature of each system, respectively. Initially, each simulated system was slowly thermalized from 10 to 300 K over 200 ps using canonical ensemble (NVT) with positional restraints of 30.0 kcal/mol·Å2 on Cα atoms of the protein. Afterward, the complex was subjected to NPT equilibration for 1300 ps with four steps of restrained MD simulations with decreasing restraints on atoms of the protein active site of 30, 20, 10 and 5 kcal/mol·Å2 and another 500 ps without any restraint. Subsequently, the entire system was run under the NPT ensemble (300 K and 1 atm) MD simulations until 200 ns was reached. The MD trajectories were collected every 10 ps. The post‑dynamic trajectories analyses in terms of RMSD, radius of gyration (Rg) and the number of atom contacts (# atom contacts) between protein and ligand were used to determine the structural variations of each complex. Additionally, the binding free energy () calculations by means of molecular mechanics/generalized Born surface area (MM/GBSA) were utilized to predict the binding affinity of the protein−ligand complexes. Whereas the per-residue decomposition free energy () was calculated to elucidate which residues of SARS-CoV-2 Mpro are important for the ligand binding. Both and were performed on 500 frames extracted from the last 20 ns of the MD production phase. It is worth nothing that the structural analysis of all MD trajectories and the binding free energies were performed using the respective CPPTRAJ utility [58] and MMPBSA.py module [59] of AMBER20.

Results and discussion

SARS-CoV-2 Mpro inhibitory activity

The twenty-one 12-dithiocarbamate-14-deoxyandrographolide analogues (3a-u) were synthesized using the standard procedure as previously reported by Arsakhant et al. [36]. The detail of general procedure for the synthesis of these analogues was described elsewhere [36]. In brief, all the analogues were synthesized and checked the purity by thin layer chromatography (TLC) and measuring the melting point. All the analogues were further identified and characterized by the nuclear magnetic resonance (NMR) spectroscopy, Fourier transform infrared (FT-IR) spectrophotometer and high-resolution mass spectroscopy (HRMS) to confirm their structures and purity. To screen the inhibitory activity of a series of C-12 dithiocarbamate andrographolide analogues toward SARS-CoV-2 Mpro, the compounds at 100 and 10 μM were tested in vitro as compared to a known SARS-CoV-2 Mpro inhibitor, rutin (positive control) [37] (Fig. 3A and B). Note that the name of each compound was also designated in the same manner as the previous study [36]. At the compound concentration of 100 μM (Fig. 3A), it demonstrated that all of the studied andrographolide analogues (3a-u) decreased the Mpro activity >70%, which was considerably larger than rutin and their parent andrographolide (1) and precursor (2). To identify the difference of inhibiory activity among tested compounds, we next performed the assay using the lower concentration of compounds at 10 μM (Fig. 3B). According to the criteria of the inhibition >50%, the result showed that compounds 3k, 3l, 3m and 3t exhibited the promising inhibitory activity against SARS-CoV-2 Mpro. However, the compound 3k was poorly soluble under the assay conditions; thus, it was not further investigated. Therefore, the remaining three compounds were then determined their IC50 and Ki values. The IC50 values of 3l, 3m and 3t were 10 ± 1, 12 ± 1 and 7 ± 1 μM, respectively, which were corresponded to the calculated Ki values of 6.8 ± 0.9, 8.3 ± 0.8 and 5 ± 1 μM (Fig. 3C). The IC50 of these compounds were also observed in the similar micromolar range with the known potent non-covalent inhibitors X77 and ML188 [60]. To gain a detailed insight on how the four compounds showed potential SAR-CoV-2 Mpro inhibition at the molecular level, the computational techniques by means of molecular docking and MD simulations were performed in the next step.

Fig. 3

Relative activity of SARS-CoV-2 Mpro in the presence of each compound at (A) 100 μM and (B) 10 μM, and (C) inhibitory profiles of compounds 3l (red), 3m (blue) and 3t (green) against SARS-CoV-2 Mpro. The Relative Fluorescence Units (RFU) produced during the initial rate period of the enzyme by time in seconds, yielding %RFU/s. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Molecular docking and rational selection of initial protein–ligand complex

To obtain more insight into how the four prominent compounds (3k, 3l, 3m and 3t) screened by in vitro assay, as shown in Fig. 3B, bind to the active site of SARS-CoV-2 Mpro, molecular docking was conducted. Initially, redocking of the crystallographic ligand back into its binding site in SARS-CoV-2 Mpro was carried out to verify the reliability of the docking simulations. The result suggested that our setting parameters of docking simulations were proper for such a system (see detail in section 2.2.2 Molecular docking of materials and methods). Afterward, each compound of interest was independently docked into the ligand-binding pocket of SARS-CoV-2 Mpro. Since the three compounds (3k, 3l and 3m) are structurally similar to each other (Fig. 2), the best pose (i.e., pose with the highest negative binding energy) of these compounds at the enzyme active site was found to be almost identical. As depicted in Fig. 4, docking results showed that the γ-crotonolactone ring of 3k, (R)-3l and 3m was well enclosed by the S1 subsite of SARS-CoV-2 Mpro formed by the side chains of F140, L141, N142 G143, S144, H163, E166 and H172. This was in good agreement with various developed peptidyl Mpro inhibitors indicating that the introduction of γ-lactam at the P1 position was strongly recognized by the S1 pocket [14], [19], [61]. Note that the binding mode of (S)-3l, (R)-3t and (S)-3t were discussed later. For the S2 subsite, i.e., the side chains of H41, C44, M49, P52 and Y54, as well as the alkyl portion of the side chains of D187 and N189, the benzyl, p‐chlorobenzyl and p-fluorobenzyl moieties of 3k, (R)-3l and 3m were deeply inserted in this hydrophobic pocket. This binding mode was also consistent with other reports demonstrating that the bulkier group such as a cyclopropyl, cyclohexyl, benzyl and fluorobenzyl at the P2 position of peptidomimetics could specifically fit with the S2 subsite, thereby improving inhibitory activities against SARS-CoV and SARS-CoV-2 Mpro [18], [19], [61]. Whereas the rest of the aryl group was favorably bound at the S4 subsite comprising the side chains of M165, P168 and N192 as well as the main chains of N189, T190 and A191. In line with our observation, several hydrophobic groups were typically introduced to interact with the S4 subsite of SARS-CoV and SARS-CoV-2 Mpro (e.g., tert-butyloxycarbonyl (BOC) and benzoxy groups), leading to the enhanced inhibitory activity of compounds against Mpro [18], [62], [63], [64]. Finally, the bulky decalin ring of all three compounds was partially occupied in the extensive S1′ subsite generated by the side chains of T25, L27, H41, C145 and H164. In the case of compound (R)-3t, it still contains the γ-crotonolactone and phenyl rings, which are expected to be a good choice to bind to the S1 and S2 subsites, respectively like in the case of the three compounds mentioned earlier. However, we found that the best pose of (R)-3t highly deviated from the presumed conformation (data not shown); thus, we alternatively selected the remaining poses that shared a similar binding mode with 3k, (R)-3l and 3m. Due to the presence of a chiral configuration at chlorcyclizine scaffold of 3l (designated as (R)-3l and (S)-3l) and 1-phenyltetrahydroisoquinoline moiety of 3t (designated as (R)-3t and (S)-3t) (see Fig. 2), the structural determinants of enantioselectivity were also explored by docking simulations. The results revealed that none of the predicted binding modes of (S)-3l could resemble the (R)-enantiomer as previously mentioned. Nonetheless, it cannot be concluded that the (S)-chlorcyclizine scaffold of 3l was not able to bind at the ligand-binding pocket in the most probable orientation as observed for (R)-3l. Therefore, to investigate the comparable effect of enantioselectivity toward protein–ligand binding, the orientation of the (S)-3l at the SARS-CoV-2 Mpro active site was prepared by manual changing in the position of chloride atom of the original (R)-enantiomer into the p-position of the benzyl ring formerly located at the S4 subsite using Discovery Studio Visualizer (BIOVIA, San Diego, CA, USA). By contrast, we found that one of the docked poses of (S)-3t reproduced a similar orientation to the selected (R)-3t, and then was chosen for comparison. In summary, our focused compounds were docked into the active site of SARS-CoV-2 Mpro and rationally selected in the most probable fashion, compared to the previously reported Mpro inhibitors, particularly at the S1, S2 and S4 subsites. These docked complexes were further applied to investigate the structural basis and dynamics behavior in aqueous solution by means of MD simulations, as provided in the following section. In addition, docking simulations of the rest of the compounds (17 compounds) were conducted to determine using a docking score alone as a cutoff value for choosing potential Mpro inhibitors or predicting the relative inhibitory activity of a compound when compared with the available experimental data. Our docking results revealed that compound 3j exhibited the highest binding affinity (AutoDock Vina score of –9.2 kcal/mol) to SARS-CoV-2 Mpro (Table S2). In contrast to in silico, this compound showed a rather low inhibitory effect on SARS-CoV-2 Mpro as compared to the more active compounds 3k, 3l, 3m and 3t mentioned above (Fig. 3B). Moreover, some of the compounds, including 3b, 3e, 3f, 3g and 3h displayed the similar docking scores (–8.5 to –8.0 kcal/mol) with our candidates 3k, 3l, 3m and 3t, although they did not show promising inhibitory activity as demonstrated by in vitro assay. Therefore, it can suggest here that only the best docking score (i.e., the most negative value in this case) for each pose may not be reliable to make a direct correlation to in vitro inhibitory activity. In fact, the reliability of using docking scores as feasible predictors for identifying SARS-CoV-2 Mpro inhibitors has been under debate from several studies [65], [66], [67]. This was the reason why more specific criteria were thus included in this study for the final selection of the docked complexes (see detail in section 2.2.2 Molecular docking). Most importantly, the evaluation of the candidates by experiment must be carried out and considered as a final decision. This would help to enhance the computational models with reliable information and give detailed insight into the activity at the molecular level.

Fig. 4

Binding conformation of each compound in the active site of SAR-CoV-2 Mpro selected from molecular docking, where the S1′, S1, S2 and S4 subsites of Mpro are represented by blue, orange, yellow and green, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Molecular dynamics of the promising compounds bound to SARS-CoV-2 Mpro

Stability of protein–ligand complexes

To assess the structural stability of each simulated model, the RMSD of complex atoms relative to the initial minimized structure was calculated and plotted along the simulation time (Fig. 5, top). RMSD profiles suggested that all complexes were likely to reach equilibrium after 120 ns with a fluctuation of ∼2.1–2.4 Å. To further support the equilibration condition of the systems and the compactness of protein structure, the Rg of protein Cα atoms was explored. The results showed that the average Rg values of all systems (Fig. 5, middle) were ∼25.4 Å, reflecting the tight compactness of the protein structure over the course of the simulation. This was also consistent with the analysis of the defined secondary structure of the protein (DSSP) framework, in which the structural features of the protein were stable along the simulation period (Fig. S3). In addition, the # atom contacts between each compound and the active site residues (Fig. 5, bottom) were calculated within 3.5 Å radius distance cutoff, i.e., any atoms closer than 3.5 Å between atoms in each compound and atoms in the active site residues were counted. The # atom contacts of these complexes appeared to fluctuate during the first 40 ns and then remained stable until the end of simulation time. The # atom contacts over the last 20 ns for 3k, (R)-3l, (S)-3l, 3m, (R)-3t and (S)-3t systems were 14 ± 4, 17 ± 4, 20 ± 4, 14 ± 4, 15 ± 4, 12 ± 4, respectively. Among all compounds, (S)-3l displayed the highest # atom contacts indicating that this enantiomer was possibly more favorable to bind to the active site residues of SARS-CoV-2 Mpro than the rest of the compounds (discussed in the following section). Altogether, based on these structural analyses, they showed that our simulation models were relatively stable; thus, the structural coordinates from the last 20 ns simulations were adopted as the production period for further analysis.

Fig. 5

Time evolution of (top) RMSD of complex atoms, (middle) Rg and (bottom) # atom contacts for 3k, (R)-3l, (S)-3l, 3m, (R)-3t and (S)-3t bound to SARS-CoV-2 Mpro during 200 ns MD simulations.

Binding affinity of protein–ligand complexes

To estimate the binding efficiency of each compound within the binding pocket of the enzyme, the ΔGbind, together with its energy components for each complex was calculated using the MM/GBSA method implemented in AMBER20. In principle, this approach combines interaction energies in the gas phase (ΔEMM) and the solvation free energy (ΔGsol) calculations. In this study, the electrostatic contribution to the solvation free energy () was computed with the generalized Born (GB) model developed by Onufriev and co-workers (igb = 2 in AMBER) [68]. Whereas the nonpolar solvation term () was assessed by the solvent-accessible surface area (SASA) with a water probe radius and a surface tension coefficient (γ) of 1.4 Å and 0.0072 kcal/(mol·Å2), respectively [69]. In general, ΔGbind of the complex is principally estimated from the difference of free energy between complex, protein and ligand. The energy components consist of ΔEMM calculated from the summation of the electrostatic (ΔEele) and van der Waals interactions (ΔEvdW), ΔGsol and entropic term (TΔS). A normal mode analysis [70] was utilized to calculate the entropy approximation by computing the translational, rotational and vibrational contributions. The free energy of the protein–ligand binding and its corresponding energy contributions of all studied complexes averaged over 500 snapshots from the last 20 ns were tabulated in Table 1. Note that due to the high computational expense of normal mode analysis, the changes in conformational entropy upon ligand binding were estimated on only 20 snapshots of the MD structures.

Table 1

MM/GBSA-based ΔGbind (kcal/mol) of each compound in complex with SARS-CoV-2 Mpro. Data are represented as mean ± standard error of the mean (SEM).

Energetics	3k	(R)-3l	(S)-3l	3m	(R)-3t	(S)-3t
Gas term
ΔEvdW	−62.52 ± 0.13	−57.91 ± 0.14	−68.57 ± 0.16	−59.87 ± 0.16	−63.50 ± 0.13	−55.40 ± 0.15
ΔEele	−14.75 ± 0.15	−17.76 ± 0.16	−22.89 ± 0.21	−7.61 ± 0.15	−16.52 ± 0.24	−15.06 ± 0.29
ΔEMM	−77.27 ± 0.20	−75.68 ± 0.21	−91.46 ± 0.29	−67.48 ± 0.20	−80.02 ± 0.30	−70.46 ± 0.34
−TΔS	26.60 ± 2.60	27.76 ± 2.66	24.30 ± 3.06	29.46 ± 3.47	25.43 ± 2.60	27.80 ± 1.86
Solvation term
ΔGsolele	44.97 ± 0.14	41.55 ± 0.14	48.56 ± 0.20	34.14 ± 0.14	45.20 ± 0.21	38.55 ± 0.22
ΔGsolnonpolar	−7.33 ± 0.02	−6.44 ± 0.01	−7.97 ± 0.02	−7.28 ± 0.02	−7.62 ± 0.02	−6.91 ± 0.02
ΔG_sol	37.64 ± 0.14	35.11 ± 0.14	40.59 ± 0.19	26.86 ± 0.13	37.58 ± 0.20	31.65 ± 0.21
Binding free energy
ΔG_bind	−13.03	−12.81	−26.57	−11.16	−17.00	−11.01

Based on the ΔEMM calculations in the gas phase, it can be seen that the main contributor to the binding of each compound to SARS-CoV-2 Mpro was the van der Waals interactions (ΔEvdW in Table 1), which were stronger than the ΔEele by ∼3- up to 8-fold. This was most likely associated with the high lipophilicity of such compound [36]. This kind of characteristic was also found in our previous MD studies and other reports of SARS-CoV-2 Mpro inhibitors, including lopinavir/ritonavir [41], peptidomimetic inhibitors (compounds N3, 11a, 13b and 14b) [71], saquinavir [72], masitinib [73] and cannabisin A [74]. By including the free energy of solvation, it was found that the nonpolar contribution () played a predominant role in the total binding free energy of all complexes, as compared to the unfavorable contribution (positive value) of the polar term (). By combing TΔS term, the calculated ΔGbind values of 3k, (R)-3l, (S)-3l, 3m, (R)-3t and (S)-3t were −13.03, −12.81, −26.57, −11.16, −17.00 and −11.01 kcal/mol, respectively. Interestingly, compound (S)-3l was predicted to have the greatest binding affinity toward SARS-CoV-2 Mpro, which was rather stronger than those of the five remaining inhibitors. On the other hand, the experimental binding free energy converted from Ki values (see section 3.1 SARS-CoV-2 M) were detected in the similar magnitude (∼−7 kcal/mol) for each system. In fact, in this study the SARS-CoV-2 Mpro inhibition assay was performed with the racemic mixture of compounds 3l and 3t, in which IC50 and Ki values would represent the overall inhibition effect from both enantiomers. Therefore, it is interesting that the compound (S)-3l should be experimentally purified and tested in vitro enzyme-based assay to validate our prediction in the future. In addition, it is worth noting that although the binding free energies based on the MM/GBSA approach did not give an absolute binding free energy in comparison with the experimental binding free energy, it is still useful to predict the trend of relative binding efficiency of the individual compound toward the target protein with fast and modest procedure. Such limitation was probably related to various factors such as GB methods, sampling protocols and simulation length [75].

Hot-spot residues upon ligand bindings

To determine the crucial residues involved in the ligand binding process, the based on the MM/GBSA method was evaluated over the same set of 500 snapshots extracted from the last 20 ns of MD simulations. The energy contributions of each residue toward the total free energy of ligand binding for each system were depicted in Fig. 6. Herein, hot-spot residues with a total energy contribution of ≤ −1.0 kcal/mol were highlighted. The binding orientations of each inhibitor in the enzyme active site were represented on the top and bottom panels of Fig. 6.

Fig. 6

The plots of calculated with MM/GBSA method for SARS-CoV-2 Mpro in complex with 3k, (R)-3l, (S)-3l, 3m, (R)-3t and (S)-3t, where the important residues involved in ligand binding highlighted in the graph (energy stabilization of ≤ −1.0 kcal/mol) were colored based on their values in the active site structures on top and bottom. The lowest to highest values were shaded from blue to red, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Among all systems, the computed suggested that the residues M49, M165 and Q189 can be considered as the important conserved residues for all ligand bindings, in which M165 showed a relatively large energy contribution ( of < –2.5 kcal/mol) for all studied complexes mainly through van der Waals interactions (Fig. 8, discussed later). This was because piperazine-dithiocarbamate moiety of 3k and (R)-3l, benzyl ring of (S)-3l, fluorobenzyl ring of 3m and 1-phenyltetrahydroisoquinoline moiety of (R,S)-3t made the hydrophobic contacts with this residue at the S4 subsite of SARS-CoV-2 Mpro. In line with this finding, the residue M165 also showed high hydrophobic interactions with peptidomimetics N3, 11a and 13b [71], boceprevir [76] and telaprevir [77]. Whereas other important residues M49 ( of < –1.8 kcal/mol) and Q189 ( of < –1.2 kcal/mol) located in the S2 subsite of the enzyme hydrophobically interacted with the decalin scaffold of all compounds as well as phenyl ring of 3k, substituted-piperazine dithiocarbamate moiety of (R,S)-3l and 3m and 1-phenyltetrahydroisoquinoline moiety of (R,S)-3t. Moreover, the additional residues (i) S46, H163 and E166, (ii) C145, H163 and R188, (iii) C44, L167, D187 and T190, (iv) S46 and D187, (v) S46, H163 and L168 and (vi) T25, H41, C44 and D187 played an essential role for the binding of 3k, (R)-3l, (S)-3l, 3m, (R)-3t and (S)-3t, respectively. These residues were also reported to be important for interacting with other SARS-CoV-2 Mpro inhibitors such as X77-mimetic candidates [78], narlaprevir [79], saquinavir, aclarubicin and faldaprevir [72].

Fig. 8

(Left) was separated into total (black bars), side chain (blue bars) and backbone (red bars) contributions for the six studied complexes of SARS-CoV-2 Mpro with (A-F) 3k, (R)-3l, (S)-3l, 3m, (R)-3t and (S)-3t, respectively. (Right) The energy contribution from electrostatic (, black line) and van der Waals (, red line) terms from the individual residue of SARS-CoV-2 Mpro. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

It should be mentioned that the phenyl ring of (R,S)-3l and (R,S)-3t and the fluorobenzyl ring of 3m were deeply inserted into the S2 subsite, as seen from the superimposition of 20 snapshots extracted every 1 ns over the last 20-ns MD simulations (Fig. 7). This was consistent well with the previous reports showing that the bulkier groups such as a cyclohexyl, cyclopropyl, phenyl ring and flurobenzyl were specific to bind to the S2 pocket of SARS-CoV and SARS-CoV-2 Mpro [80]. In addition, γ-crotonolactone ring of most compounds somewhat fitted well into the S1 subsite, as the chemical structure of this moiety was highly similar to the cyclized glutamine mimetic (e.g., (S)-γ-lactam group), which was found to specifically bound to the S1 pocket of SARS-CoV-2 Mpro [19]. Nevertheless, the γ-crotonolactone ring of (S)-3t slightly moved outward the S1 subsite, resulting in the less tight packing within this pocket and probably leading to the lower binding affinity of (S)-3t than that of (R)-3t (Table 1). For the S4 subsite, the benzyl, chlorobenzyl, fluorobenzyl and tetrahydroisoquinolin moieties of the respective compound 3k, (S)-3l, 3m and (R,S)-3t were surrounded by this shallow hydrophobic subsite. In correspondence with several X-ray structures, they demonstrated that the S4 subsite of SARS-CoV and SARS-CoV-2 Mpro was more favorable to bind to various bulkier and hydrophobic groups such as benzoxy [63], BOC [18], BOC-serine [81] and indole groups [19]. Instead, the chlorobenzyl of (R)-3l was more likely to bind to the S2 subsite rather than the S4 subsite, as previously observed in other cases mentioned earlier. This might be the reason why (R)-3l showed much weaker binding efficiency as compared to its S configuration (Table 1). Unexpectedly, it can be seen from the overlay structures of each compound that the decalin ring of all compounds could not fit well within the S1′ subsite or even move out from this subsite of SARS-CoV-2 Mpro, when compared with the starting structure taken from the molecular docking (Fig. 4). Indeed, the ligand conformations taken from the MD simulation were likely to drift from their initial conformation to their relaxed structures and avoid the bad contacts between the active site residues and ligand. In such cases, the modification of decalin moiety with O-acetylation on the two hydroxyl groups could lead to the steric hindrance with the S1′ subsite residues and move outward to the solvent-expose area in our MD simulations. The obtained information was supported by a high fluctuation of the decalin moiety detected in all systems over the last 20-ns simulation (Fig. 7). Therefore, we suggest modifying this moiety back to the original one to reduce such steric effect together with decreased molecular weight and improved drug-likeness [82], and may also enhance the hydrogen bond (H-bond) formation with the residues at the S1′ subsite.

Fig. 7

Overlay structures over the 20 snapshots of each focused compound in complex with SARS-CoV-2 Mpro derived from the last 20-ns MD simulations, where the S1′, S1, S2 and S4 subsites of Mpro are shaded by blue, orange, yellow and green, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

To further evaluate the energetic contributions, the degree of stabilization from individual residue marked in Fig. 6 was considered the energy contributions from its backbone and side chain atoms, as shown in Fig. 8(left). Additionally, the contributed energies from electrostatic () and van der Waals () terms were given in Fig. 8(right). The results suggested that most of the critical residues for each complex had a tendency to stabilize all compounds through their side chains, as indicated by the higher energy contribution (larger negative values) from the side chain atoms (blue bars in Fig. 8, left) than those of the backbone atoms (red bars). We also found that the main energy contribution for binding of compounds with each residue came from the van der Waals interactions. It can be seen from the greater negative values (red lines in Fig. 8, right) up to ∼−5.5 kcal/mol, as opposed to the electrostatic interactions (∼−0.6 to ∼4.0 kcal/mol, black lines), which were in good accordance with the ΔEMM results (Table 1) as well as the very low H-bond occupations detected in all complexes (Table 2 and Fig. 9, discussed later). The van der Waals contribution from the three overlapping residues M49, M165 and Q189, pivotal for binding of all compounds, was lower than −2 kcal/mol in all systems, reflecting the importance of these residues for recognizing the ligand bindings. This finding was consistent well with our previous MD studies on the HIV-1 protease inhibitors lopinavir and ritonavir [41] and peptidyl inhibitor N3 [71] against SARS-CoV-2 Mpro. It was also supported by the computational alanine scanning calculation suggesting that mutations of these residues into alanine culminated in the weakened binding interactions with masitinib and α-ketoamide inhibitor 13b [73].

Table 2

Percentage of H-bond occupations detected between the heteroatoms of each compound and the residues of SARS-CoV-2 Mpro calculated over the last 20 ns of MD simulations. Note that hydrogen bonding >20% was reported here.

H-bonding	H-bond occupations (%)
	3k	(R)-3l	(S)-3l	3m	(R)-3t	(S)-3t
O6⋯H − OG1(T24)	–	–	28.3	–	48.9	–
O6⋯H − N(T24)	–	–	27.4	–	–	57.4
O6⋯H − OG1(T25)	–	–	25.0	–	–	82.8
F1⋯H − ND1(H41)	–	–	–	50.4	–	–
O1⋯H − NE2(H163)	–	75.9	–	–	–	–
O2⋯H − NE2(H163)	–	48.2	–	–	20.5	–
O1⋯H − N(E166)	–	–	95.0	29.0	–	–
O2⋯H − N(E166)	–	–	–	54.5	–	–

Fig. 9

Representative 3D structures showing H-bond formations (black dashed lines) of each focused compound with the active site residues of SARS-CoV-2 Mpro.

Intermolecular hydrogen bonds between protein and ligand

As stated in Table 1, the main contribution for such andrographolide analogues binding to SARS-CoV-2 Mpro was the van der Waals interactions. Thus, a low number of H-bond formations between each compound and the binding residues were expected. To measure such interaction, the percentage of H-bond occupations between each compound and the SARS-CoV-2 Mpro residues during the last 20-ns simulations were determined based on the two following geometric criteria: (i) an acceptor⋯donor distance of ≤ 3.5 Å and (ii) acceptor⋯H-donor angle of ≥ 120°. The percentage of H-bond occupations was given in Table 2, whereas the H-bond patterns of the compounds in the binding pocket of SARS-CoV-2 Mpro were illustrated in Fig. 9. Expectedly, most compounds did not form strong H-bonds (>80% occupation) with the active site residues. This was except for the E166 backbone nitrogen, which formed H-bond with the O1 atom (see Fig. 2 for atomic labels) of (S)-3l at 95.0%. It has been reported from several investigations that this residue could interact with the peptidomimetic inhibitors (N3, 11a and 13b) [71], nirmatrelvir [83] and boceprevir [76] via H-bond interactions. We also found that the hydroxyl group of T25 side chain formed a strong H-bond with the O6 atom of (S)-3t at 82.8%. However, no H-bond interactions between 3k and its surrounding residues were detected. Whereas the remaining compounds formed weak (<50% occupation) and/or moderate (50–79% occupation) H-bond interactions with their surrounding residues at the enzyme active site (Table 2). Interestingly, it should be emphasized that the residues H163 and E166 also partially participated in stabilizing γ-crotonolactone ring of (R,S)-3l, 3m and (R)-3t, which was similar to (S)-γ-lactam ring of the previously reported peptidyl aldehyde and α-ketoamide inhibitors [18], [19], enclosed by the S1 subsite of SARS-CoV-2 Mpro. Taken altogether, it can be concluded that H-bond formations played a minor role in binding recognition between SARS-CoV-2 Mpro and andrographolide analogues. By contrast, the van der Waals interactions contributed from each binding site residue of SARS-CoV-2 Mpro were considered as the main contributor to the binding of all compounds, as shown in red lines in Fig. 8(right).

Solvent accessibility at the SARS-CoV-2 Mpro active site

The effect of water accessibility on the SARS-CoV-2 Mpro active site upon ligand binding was measured by calculating SASAs on the amino acid residues within 5 Å of each compound (Fig. 10A). Note that only protomer A contained the ligand inside the enzyme active site. The time evolution of SASA for each studied system was plotted in Fig. 10B, and the average SASA values derived from protomer A and protomer B from the last 20 ns of MD simulations were also listed in the green and yellow texts, respectively. The obtained results showed that the average SASAs for the apo form (without ligand bound, protomer B in Fig. 10A) of 3k, (R)-3l, (S)-3l, 3m, (R)-3t and (S)-3t system were 1052.00 ± 99.58, 1029.06 ± 64.88, 1053.62 ± 58.76, 1054.01 ± 81.74, 1076.63 ± 48.55 and 950.19 ± 45.16 Å2, respectively. Upon ligand binding in the enzyme active site, i.e., protomer A, the average SASAs of all simulated models diminished by ∼120 to 450 Å2. This event was in accordance with our previous MD reports on lopinavir/ritonavir [41], peptidomimetic inhibitors (compound N3, 11a, 13b and 14b) [71] and halogenated baicalein [84] bound to SARS-CoV-2 Mpro, which showed a significant reduction in SASAs during the binding process. Moreover, to compare the enantiomer system, it can be noticed that the (S)-3l system exhibited lower SASAs than the (R)-3l system, whereas the (S)-3t system showed higher SASAs than the (R)-3t system. This finding implied that (S)-3l and (R)-3t could bind to the binding pocket of SARS-CoV-2 Mpro better than their opposite chiral configuration. This was supported by the ΔGbind calculations (Table 1), indicating the stronger binding affinity of (S)-3l than (R)-3l, and (R)-3t than (S)-3t toward SARS-CoV-2 Mpro.

Fig. 10

(A) SARS-CoV-2 Mpro homodimer in the presence and absence of ligand in protomer A (green) and protomer B (yellow), respectively. Note that the selected residues within 5-Å sphere from the ligand were utilized for SASA calculations. (B) SASA plots of the six studied systems during MD simulations. The orange and green texts represented the average SASA values (mean ± SD) for the protomer A and protomer B derived from the last 20 ns of individual MD simulations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Conclusion

In this study, the combined experiment and computational techniques were conducted to discover a novel SARS-CoV-2 Mpro inhibitor based on our recent series of 12-dithiocarbamate-14-deoxyandrographolide analogues. Twenty-one andrographolide analogues (3a-u), including their parent compound (1) and precursor (2), were tested in vitro inhibitory activity against SARS-CoV-2 Mpro. The results revealed that compounds 3k, 3l, 3m and 3t showed promising inhibitory activity toward SARS-CoV-2 Mpro better than the known inhibitor rutin. To further clarify the mechanism of action of these compounds at the molecular level, molecular docking and MD simulations were carried out. Initially, the binding orientations of each compound including the enantiomer of 3l and 3t ((R,S)-3l and (R,S)-3t), bound to the active site of SARS-CoV-2 Mpro were studied by molecular docking. The optimum structures for each complex were selected based on the similar orientations of each compound moiety to the preferential scaffold of several reported SARS-CoV-2 Mpro inhibitors, i.e., γ-crotonolactone ring and phenyl ring should bind to the S1 and S2 subsites, respectively. Afterward, 200-ns MD simulations were performed to investigate the stability of the docked complexes. The structural analyses in terms of RMSD, Rg and # atom contacts indicated that our complex models were stable along the simulation time. The MM/GBSA-based ΔGbind calculations suggested that the van der Waals interactions were the primary driving force for the molecular complexation between all studied compounds and SARS-CoV-2 Mpro. Moreover, it was found that compound (S)-3l exhibited the highest binding affinity toward SARS-CoV-2 Mpro than those of the remaining compounds. The hot-spot residues, including T25, H41, C44, S46, C145, H163, E166, L167, D187, R188 and T190, and particularly for the three residues M49, M165 and Q189 identified by calculations were considered as the key residues responsible for ligand binding. However, it can be noticed from the overlay structures over the last 20 ns that the acetoxy groups on the decalin rings of each analogue may cause the steric hindrance with the residues at the S1′ subsite and move out from this pocket to the solvent-expose area. This finding could help to guide the synthesis of a new series of andrographolide analogues with enhanced inhibitory activity in the future. Taken together, the experimental and theoretical results presented here provide important information to describe the structure-activity relationship and the possible binding interactions of the four prominent andrographolide analogues to SARS-CoV-2 Mpro, which could also be useful for further optimizations of more potent anti-SARS-CoV-2 Mpro.

CRediT authorship contribution statement

Bodee Nutho: Conceptualization, Investigation, Data analysis, Writing - original draft, Writing - review & editing, Funding acquisition, Project administration. Patcharin Wilasluck: Investigation, Data analysis. Peerapon Deetanya: Investigation, Data analysis. Kittikhun Wangkanont: Investigation, Data analysis, Resources, Writing - review & editing. Patcharee Arsakhant: Investigation, Data analysis. Rungnapha Saeeng: Investigation, Data analysis, Resources, Writing - review & editing. Thanyada Rungrotmongkol: Conceptualization, Resources, Software, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.