A molecular modelling approach for identifying antiviral selenium-containing heterocyclic compounds that inhibit the main protease of SARS-CoV-2: an in silico investigation

Abstract

Coronavirus disease 2019 (COVID-19), an infectious disease caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has been declared a global pandemic by the World Health Organization, and the situation worsens daily, associated with acute increases in case fatality rates. The main protease (M^pro) enzyme produced by SARS-CoV-2 was recently demonstrated to be responsible for not only viral reproduction but also impeding host immune responses. The element selenium (Se) plays a vital role in immune functions, both directly and indirectly. Thus, we hypothesised that Se-containing heterocyclic compounds might curb the activity of SARS-CoV-2 M^pro. We performed a molecular docking analysis and found that several of the selected selenocompounds showed potential binding affinities for SARS-CoV-2 M^pro, especially ethaselen (49), which exhibited a docking score of −6.7 kcal/mol compared with the −6.5 kcal/mol score for GC376 (positive control). Drug-likeness calculations suggested that these compounds are biologically active and possess the characteristics of ideal drug candidates. Based on the binding affinity and drug-likeness results, we selected the 16 most effective selenocompounds as potential anti-COVID-19 drug candidates. We also validated the structural integrity and stability of the drug candidate through molecular dynamics simulation. Using further in vitro and in vivo experiments, we believe that the targeted compound identified in this study (ethaselen) could pave the way for the development of prospective drugs to combat SARS-CoV-2 infections and trigger specific host immune responses.

Open in new tabDownload slide

INTRODUCTION

In December 2019, a cluster of pneumonia cases of unknown aetiology was identified in Wuhan, a city in the Hubei province of China. The patient’s symptoms were largely mild-to-moderate and included fever, dry cough and difficulty breathing, with an excessive amount of mucus. Cardiovascular and neurological complications were observed in severe cases (1–3). This acute respiratory tract infection was caused by a novel coronavirus, called severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2, formerly 2019-nCoV). The World Health Organization (WHO) named the disease coronavirus disease 2019 (COVID-19) and declared the outbreak a Public Health Emergency of International Concern (PHEIC) on 30 January 2020. The symptoms identified were similar to those of previous human disease outbreaks caused by coronaviruses, including severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV). However, the mortality rate of SARS-CoV-2 has been reported at 3%, whereas the mortality rates for SARS-CoV and MERS-CoV were reported to be 9.6% and 34%, respectively (4). On 11 March 2020, the WHO declared the COVID-19 outbreak a global pandemic and, as of 10 October 2020, more than 36.9 million cases have been reported, resulting in over 1.07 million deaths (5,6).

SARS-CoV-2 has been identified as a member of β-coronavirus clade b, and its genome consists of a non-segmented, single-stranded RNA. Generally, a large polyprotein, 1ab (replicase 1ab) with a molecular weight of 800 kDa, is produced by the gene transcription of β-coronaviruses, and the proteolytic cleavage of this polyprotein produces several proteins.

Recently, Liu and his research group confirmed the presence of the main protease (M^pro) enzyme in SARS-CoV-2 (7). Protein synthesis is dependent on the activity of M^pro, also known as 3CL^pro, which, along with papainlike proteases, regulates polyprotein processing. M^pro activity is regarded as a prerequisite for viral replication; hence, viral replication would be blocked by inhibiting the activity of M^pro [7]. Interestingly, SARS-CoV-2 proteins have been shown to interfere with host immune responses, and M^pro-specific T cells have been identified in patients infected with SARS-CoV-2 (8,9). Additionally, previous studies have demonstrated that the papainlike proteases produced by SARS-CoV and MERS-CoV can hinder the immune response, and SARS-CoV antagonises the antiviral response mediated by interferons (10). Moreover, studies have previously reported that the reproduction rate of SARS-CoV-2 is significantly higher than that of its predecessors (11). An initial estimate by the WHO suggested a basic reproduction number (R₀) of 1.4–2.5 (average 1.95); however, a much higher R₀ value of approximately 5.7 (with a 95% confidence interval of 3.8–8.9) was calculated in a later study (12). We hypothesise that SARS-CoV-2 M^pro may be responsible for curbing immune functions and that the inhibition of SARS-CoV-2 M^pro activity might increase immune function, thereby increasing immunity.

Currently, no established therapies exist for the treatment of COVID-19, and researchers continue to search for suitable therapeutic agents for acute disease. The treatment measures that are currently implemented are insufficient to avoid the risk of target organ damage. Recently, many clinical trials have investigated various antiviral agents and active chemical entities that have been repurposed from the treatment of other viral infections. However, the results have thus far been inconclusive, and the findings are too preliminary for use in the development of pharmaceutical products outside of the clinical trial setting (13,14). Therefore, further investigation remains necessary to support the development of pharmacologically active substances to combat COVID-19.

Ample mechanistic clinical data have shown that metals and their compounds play important roles in immune functions. Trace elements, including iron (Fe), zinc (Zn), manganese (Mn) and copper (Cu), regulate various important and complementary functions, such as DNA replication, and these metals operate as co-factors for several enzymes and co-enzymes. Recently, the attachment of cobalt-III was reported to inhibit the active site of SARS-CoV-2 M^pro (15). Importantly, the biological roles played by selenium (Se) were revealed when glutathione peroxidase (GPx) was found to contain Se in its active site. Se is an essential element in the antioxidant defence system, in various metabolic pathways and, most importantly, in the function of both the innate and adaptive immune systems (16).

Se plays diverse roles in the regulation of immune function and is detectable in various tissues, including lymph nodes, liver and spleen. In addition, multiple cellular and viral mechanisms involving Se and selenoproteins, such as thioredoxin reductases, have been shown to indirectly regulate the activation of the immune response (17). Increased Se concentrations enhance cytosolic and phospholipid GPx and decrease viral pathogenicity, and several viruses have been reported to be more virulent in selenium-deficient hosts because Se is involved in host immune responses (18,19). In mammalian redox biology, Se is a necessary trace element, and a deficiency in dietary Se can cause oxidative stress in the host, resulting in the alteration of the viral genome, which can trigger a normally benign pathogenic virus to become highly virulent (19,20). Se is involved in the protection against many evolutionarily distinct viruses via potential immunomodulatory effects, which is consistent with a role for Se in the immune system (21).

Se deficiency can trigger host immune system impairment and increase viral pathogenesis by causing alterations in the viral genome, resulting in the rapid and reproducible introduction of genetic mutations in benign RNA viral variants (influenza, poliovirus and coxsackievirus) that transform them into virulent variants (22–24). In chickens, Se can induce an immune response to a live bivalent infectious bronchitis coronavirus vaccine due to synergistic effects with the saponins of ginseng stem leaves (25). These clues suggest that this trace element may be an effective component in the treatment for SARS-CoV-2 infection. In the present study, we assessed the ability of several Se-containing heterocyclic compounds to inhibit SARS-CoV-2 M^pro and demonstrated the biological activities of these compounds using computational biology techniques.

MATERIALS AND METHODS

Retrieval of the protein sequence and homology modelling

The sequence of the SARS-CoV-2 M^pro protein was obtained from the National Center for Biotechnology Information (NCBI) database (PDB: 6M03_A, GI: 1820435677), as the query sequence. The tertiary structure of SARS-CoV-2 M^pro was modelled based on the available template deposited in the Protein Data Bank (PDB, https://www.rcsb.org/). Nine templates for SARS-CoV-2 M^pro, with high resolution (PDB ID: 6WQF, 5R80, 5RE6, 5REA, 6 LU7, 6 M03, 6 W63, 6XCH and 6WTT), were downloaded from the PDB database. Homology modelling was performed using EasyModeller, a graphical user interface version of the program Modeller (26,27). EasyModeller is completely automated for the generation of energy-minimised protein models, and ModRefiner was used to minimise and correct the structure (28).

Model evaluation and validation

The stereochemical accuracy of the model’s protein structure and its overall structural geometry were confirmed using the Procheck program (29). Ramachandran plot statistics were analysed to evaluate the stability of the model and to confirm the residues. The overall quality of the model was measured using the ProsaWeb server (30). We also performed Ramachandran plot and Z-score analyses for four selected high-resolution SARS-CoV-2 M^pro structures obtained from the PDB database, including 6 LU7 (Resolution: 2.16 Å), 6 W63 (Resolution: 2.1 Å), 6WQF (Resolution: 2.3 Å) and 6XCH (Resolution: 2.20 Å). In addition, to understand the stable nature and flexibility of the SARS-CoV-2 M^pro model, normal mode analysis (NMA) was performed (31,32). The iMod tool was utilised to determine the stability level; this tool can also function as a valid alternative to atomistic simulation (33–35). The motion stiffness was presented through a given value, and the elastic network model, deformability, eigenvalue and covariance matrix were calculated.

Molecular docking analysis

Ligand retrieval and preparation

A total of 49 Se-containing heterocyclic compounds were selected for the present study (36). Compounds 1–42 were drawn using ChemDraw version 16.0 (PerkinElmer ChemDraw Professional) and compounds 43–49 were obtained from the PubChem database (CID No.: 9990475, 24786497, 11539212, 134741, 128911, 10544652, 10387485; Supplementary Table S1). The 3D structures of the ligands were built using PyRx 0.8 (37). The energy form of all compounds was minimized and converted into pdbqt format by Open Babel in PyRx version 0.8.

Protein preparation

We used the protein preparation wizard of UCSF Chimera version 1.11.2 to prepare the modelled SARS-CoV-2 M^pro protein for docking analysis by adding hydrogens and Gasteiger–Marsili charges (38,39). The prepared file was then converted into pdbqt format using Open Babel (40).

Molecular docking

Docking calculations were performed using AutoDock Vina in PyRx, version 0.8. The parameters used for the docking simulation were set to default. During the course of the docking procedure, the grid box was centred to cover the binding site residues and permit the ligand to move freely. The box was set to 32.00499 × 27.03345 × 29.09178 Å. In addition, we also docked the compound with the highest docking score with several proteins from the SARS-CoV-2 proteome, including, the nucleocapsid (N) protein, the receptor-binding domain (RBD) of the spike (S) protein, papainlike protease (PL^pro), RNA-dependent RNA polymerase (RdRp), nsp9 RNA-replicase, SARS-CoV-2 helicase (nsp13) and uridylate-specific endoribonuclease (nsp15). For docking analysis, the following receptors were selected: PDB ID 6M3M for N protein, PDB ID 2GHV for the RBD of the S protein, PDB ID 6W9C for PL^pro, PDB ID 6 M71 for RdRp, PDB ID 6ZSL for SARS-CoV-2 helicase (nsp13), 6WC1 for nsp9 RNA-replicase, and 6XDH for nsp15. The amino acid residues were kept flexible. AutoDock Vina was implemented via the shell script offered by the program’s developers. Docking results were reported as a negative score, in kcal/mol, with the lowest docking score indicating the highest binding affinity (41,42).

Validation of docking approach

We repeated the molecular docking analysis for validation of the strategy of the docking procedures for the best ligand and the respective SARS-CoV-2 proteins. For validation purposes, first we used GC376 as positive control and a recent experiment reported that GC376 bound with the active site of the SARS-CoV-2 M^pro (43). Afterwards, we have performed cross-docking analysis with selected compounds with one of the crystal structure of SARS-CoV-2 M^pro (PDB ID: 6 LU7) available in PDB database. In addition, as no respective ligand molecules are available in the PDB database, we selected zidovudine (PubChem CID: 35370), coenzyme-A (PubChem CID: 87642), elbasvir (PubChem CID: 71661251), remdesivir (PubChem CID: 121304016) as positive controls for N protein, RBD, PL^pro and RdRp, ioxilan for nsp9, cangrelor for nsp13, and citric acid for nsp15, respectively (44–49).

Ligand validation

The selected selenocompounds were subjected to biological activity calculations using an online validation tool, the Molinspiration Cheminformatics server (www.molinspiration.com) (50). Molinspiration calculates the molecular properties associated with drug-likeness and bioactivity, including predictions regarding whether a compound is a nuclear receptor ligand, a G protein-coupled receptor (GPCR) ligand or an enzyme inhibitor. The calculated drug-likeness scores of all compounds were compared with those of the positive control GC376.

Ligand-based ADME prediction

The evaluation of the pharmacokinetic properties of all identified compounds was carried out using Lipinski’s rule of five (51). According to this rule, a compound could show optimal druglike behaviour if it fulfils at least four of the five characteristics, namely molecular weight < 500 Daltons, ≤5 hydrogen bond donors, ≤10 hydrogen bond acceptors, lipophilicity <5 and molar refractivity between 40 and 130. The web tool SwissADME, a convenient tool in drug discovery, was used to analyse these properties for all of the identified compounds (52). Compounds that meet Lipinski’s condition are considered ideal drug candidates.

Prediction of toxicological properties

The toxicological properties of the selected selenocompounds were also predicted because toxicity is a major concern for the administration of any drug. The admetSAR server (http://lmmd.ecust.edu.cn/admetsar1/predict/), which was previously described as an important and convenient prediction tool for drug discovery, was used for this purpose (53). Ames toxicity, carcinogenicity, biodegradation, acute oral toxicity, rat acute toxicity and inhibitory effects on human ether-a-go-go-related gene (hERG) were predicted.

Molecular dynamics simulation

The molecular dynamics simulation study was implemented in the YASARA dynamics commercial software package to validate the docking study (54). The AMBER14 force field was used, with an initial energy minimisation process, and was implemented with the steepest gradient approach (5000 cycles) (55). The system was neutralised with the addition of water molecules and 0.9% NaCl at 310 K. The particle mesh Ewald method (PME) was applied to calculate long-range electrostatic interactions (56). A cubic simulation cell was created, which was 20 Å larger than the complexes, in all cases. A periodic boundary condition was maintained, and the temperature of the system was controlled by a Berendsen thermostat. The time step of the simulation system was set to 1.25 fs (57). Finally, the simulation was performed for 100 ns, and after every 100 ps, the simulation trajectories were saved. The simulation trajectories were analysed to calculate root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA) and hydrogen bonds (58–60).

To calculate the binding-free energy, the Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) method was applied (61). The YASARA macrofile was edited and developed to calculate the binding-free energy. The following equation was employed to estimate the free energy:

ΔGbind = ΔGcomplex_(minimised) – [ΔGligand_(minimised) + ΔGreceptor_(minimised)].

ΔGbind = ΔG_MM + ΔG_PB + ΔG_SA − TΔS.

Here, ΔG_MM is the sum of van der Waals and electrostatic interactions; ΔG_PB and ΔG_SA are the polar and non-polar solvation energies, respectively; and TΔS is the entropic contribution (62).

In silico IC50 calculation

Executing AutoDock

The inhibition constant in this study was determined by AutoDock 4.0 tools. AutoDock software calculates and predicts the interactions between a ligand molecule and a protein molecule, based on predefined parameters (63). The interactions between the molecules can be calculated at a specified user-defined region of the protein. This region must be defined by the user through the Grip map option or GridBox. Thus, the use of the GridBox at the binding site or active site, or at other essential protein regions, is essential to perform AutoDock analysis. Before executing the AutoDock, the ‘.pdb’ files of the protein and ligand have to be moved into one folder.

Analysis in AutoDock can be divided into the following categories: (a) initialising molecules and (b) running AutoGrid.

Initialising molecules

Initialising the molecule involved the addition of hydrogen atoms. The ligand molecule required the addition of the Gasteiger charge, the identification of aromatic carbons, the detection of rotatable bonds and the setting of the torsional degrees of freedom (TORSDOF) value. The protein must be initialised manually, whereas the ligand is automatically initialised when opened in the tool. The receptor and ligand files were saved in .pdbqt format. Next, the receptor was opened again using the ‘Grid’ menu and the ‘Macromolecules’ sub-menu, and ‘open’ was selected. The ligand opens using the ‘Set Map Types’ sub-menu and ‘Open Ligand’. Then, the GridBox is set in AutoDock to cover the identified binding sites. AutoDock only analyses the ligand molecule’s interactions and the amino acids that are present within the GridBox. The GridBox’s size can be increased or decreased using the number of points in the X/Y/Z dimension. The position of the GridBox can be adjusted to cover the binding site or binding residues, using the ‘Center GridBox’ field that moves the GridBox in the X-, Y- and Z-axis.

Running AutoGrid and AutoDock

After setting up the GridBox, the AutoGrid file was saved. When in the ‘Grid’ menu, select the ‘Output’ and ‘Save GPF’ sub-menu. The file saves in .gpf format. Meanwhile, the AutoDock file saving procedure is to choose the ‘Docking’ menu and select ‘Lamarckian GA’ in the ‘Output’ menu. The file saves in .dpf format. The AutoGrid and AutoDock programs were run through a command prompt directed to the folder to be docked using the command: autogrid4 -p control.gpf -l control.glg. After the AutoGrid calculation was completed successfully, AutoDock was run using the command: autodock4 -p control.dpf -l control.dlg. The results of the docking calculations are obtained in a notepad file format, reporting the values for binding energy and the estimated inhibition constant (64).

RESULTS

Tertiary structure of SARS-CoV-2 M^pro

The model was based on the template sequences described in the Materials and Methods section. Homology modelling was performed using EasyModeller, which utilises the building modules of Modeller to generate 3D models. Three models of SARS-CoV-2 M^pro were generated, all of which had identical structural characteristics to the query protein. The best model [the model with the lowest discrete optimised protein energy (DOPE) score] was selected for further analysis. The structure of the model was refined and used for further analysis (Figure 1A).

Quality of the predicted model

Different structure validation programs are used for the evaluation of the generated model proteins, which included geometrical conformations and stereochemical quality assessments. The PROCHECK program was used to perform the Ramachandran plot calculations. The Ramachandran plot of the model protein indicated that 95.1% of residues were in the most favourable region, 3.4% were in the allowed region and 0.4% were in the disallowed region (Figure 1C). In addition, the ProsaWeb-derived Z-score of the model was −7.32, indicating that the model represents a nuclear magnetic resonance (NMR)-quality structure (Figure 1B). In addition, the Rama favoured region for the crystal structures of SARS-CoV-2 M^pro were shown to be 90.6% for 6 LU7, 91.7% for 6 W63, 87.9% for 6WQF and 88.3% for 6XCH (Supplementary Figure S1). 6 LU7 showed a Z-score of −7.34, which was the highest among the selected crystal structures, including the model SARS-CoV-2 M^pro. The stability and the quality of the hypothetical SARS-CoV-2 model were assessed through B-factor, eigenvalue, covariance map, deformability and elastic network model. Figure 2 revealed an insignificant hinge and an average RMS in the B-factor. A low chance of deformation for the predicted homology model was indicated by a high eigenvalue, 1.382697e⁻⁴, and the correlation and elasticity demonstrated the high quality of the modelled SARS-CoV-2 M^pro.

Computational molecular docking studies

Molecular docking was performed using AutoDock Vina in PyRx, version 0.8, to verify the binding model of the selenocompounds with our hypothesised protein. Along with GC376, most of the selected selenocompounds interacted with the SARS-CoV-2 M^pro protein model. However, only 49 showed a higher affinity toward the receptor than the control, GC376. The docking score for 49 was −6.7 kcal/mol, whereas the control compound, GC376, possessed a docking score of −6.5 kcal/mol (Table 1). In addition, the visualisation of the docking experiment findings for 49 revealed hydrophobic interactions with Thr25, Leu27 and Met49 residues and attractive charges attributed to the Met49 and Cys145 residues (Figure 3A). In addition, visualisation revealed that the positive control, GC376, also yielded a hydrophobic (pi–alkyl) interaction with the Met49 residue (Figure 3B). The scores for both 26 and 6 were − 6.4 kcal/mol, predicting a high binding affinity. In this experiment, we selected the 16 selenocompounds with the best docking scores for further analysis.

The 16 selective compounds were further docked against the crystal structure of SARS-CoV-2 M^pro (PDB ID: 6 LU7) (Table 2). Compounds 49, 26 and 44 showed high affinities towards the receptor, with docking scores of −7.7, −7.2 and − 7.4 kcal/mol, respectively, compared with a docking score of −8.1 kcal/mol for GC376. The visualisation of 26 depicted the formation of a hydrogen bond with Met49, hydrophobic interactions with His41 and Cys145, and attractive charges attributed to Met49 and Cys145, with the absence of unfavourable bumps. Compound 44 formed a hydrogen bond with Arg188, hydrophobic interactions with Met49, Thr190, Glu166, Pro168 and Gln189, and attractive charges were attributed to Met165, with the absence of unfavourable bumps. In the case of 49, a hydrogen bond was formed with Ser144 and Cys145, and hydrophobic interactions were observed for His41, Leu141, Asn142, Gly143, Cys145, Met165 and Glu166 (Figure 4A). GC376 formed hydrogen bonds with Phe140, Met165, Glu166 and Arg188 and had hydrophobic (pi–alkyl) interactions with His41 and Pro168 residues. GC376 formed unfavourable interactions with His163, Thr190 and Gln192, and attractive charges were attributed to Glu166 (Figure 4B).

In addition, the compound with the best docking score against both the model and crystal structures of SARS-CoV-2 M^pro, 49, also interacted with the N protein, the RBD of the S protein, nsp9 RNA-replicase, SARS-CoV-2 helicase (nsp13) and uridylate-specific endoribonuclease (nsp15), with docking scores of −8.4, −7.8, −6.7, −7.6 and − 6.9 kcal/mol, respectively (Table 3). Intriguingly, the docking scores were higher than those for the controls that were used for the respective receptor proteins: as zidovudine showed a docking score of −7.1 kcal/mol when interacting with the N protein; co-enzyme A exhibited a docking score of −6.4 kcal/mol during the interaction with the RBD of the S protein; ioxilan showed a docking score of −4.3 kcal/mol when interacting with nsp9; cangrelor had a docking score of −7.5 kcal/mol when interacting with nsp13 and citric acid exhibited a docking score of −5.5 kcal/mol during the interaction with nsp15 (Table 3). However, 49 showed less binding affinity towards PL^pro and RdRp than the controls (Table 3). The results of the docking experiments are shown in Figures 3–11, Supplementary Figures S1–S7, Tables 1–3 and Supplementary Table S2.

Biological activities of the drug candidates

The selected compounds were screened in parallel with GC376 using the Molinspiration Cheminformatics server to predict their biological activities. The findings demonstrated that 8 possessed greater enzyme inhibitor activity than GC376. In addition, 44 exhibited the highest nuclear receptor ligand activity among the selected compounds, which was higher than that for the standard GC376. Compounds 49, 44 and 26 displayed the highest kinase-inhibiting activity among the targeted compounds, which were higher than that of GC376; however, GC376 showed the best protease inhibitor activity. The results are shown in Table 4.

Drug-likeness of the potential selenocompounds

SwissADME, an online-based tool, was used to determine the drug-likeness properties of the selected selenocompounds with the highest binding affinities. Interestingly, all of the targeted compounds followed both Lipinski’s rule of five and Veber’s rule. The ADME (absorption, distribution, metabolism and excretion) properties are shown in Table 5.

Toxicological properties of the potential selenocompounds

The toxicological characteristics of the 16 selected compounds were predicted using the online toolkit admetSAR (Table 6). The results indicated that all of the selenocompounds were likely non-carcinogenic, and none were found to be biodegradable. Importantly, the Ames toxicity test was positive for only 13 and 26, indicating that most of the compounds posed no threat of mutagenicity. Additionally, all of the selected selenocompounds were found to be weak hERG inhibitors. All of the compounds exhibited weak rat acute toxicity, with a median lethal dose (LD₅₀) ranging from 2.1398 to 2.5469 mol/kg. The acute oral toxicity values for all of the selected compounds lie within Category III, which includes compounds with LD₅₀ values greater than 500 mg/kg but less than 5000 mg/kg and are usually considered suitable as drugs (65). Thus, all of the examined selenocompounds, except 13 and 26, passed all of the selected toxicological parameters.

Molecular dynamics simulation

The RMSD for both the control and drug–protein complexes was evaluated to understand their stability levels. We used 49 and the model SARS-CoV-2 M^pro complex as the inputs. Figure 12A shows that both the control and drug complex exhibited an upper trend from the starting point of the simulation; however, the drug complex had increased fluctuations compared with the control at the 20–30 ns period. Thereafter, the control and drug complex tended to stabilise and exhibited a steady state, denoting the complex’s integrity over the simulation trajectory. Moreover, the solvent-accessible surface areas (SASAs) of both complexes were analysed to understand the changes in the protein volume. Interestingly, the control complex featured an expansion of the surface area, with a higher SASA value than the drug complex (Figure 12B). The drug complex remained in a stable state, despite some fluctuation points. The result from the SASA descriptors demonstrated that the drug complex did not increase in volume during the simulation process, whereas the control drug appeared to be more unstable.

The degree of protein firmness correlates with the Rg, where a higher Rg profile indicates the loose packaging system, and a lower Rg profile is associated with a more compacted protein formation. Figure 12C indicates that the drug complex had a higher Rg value than the control, which might be responsible for the more mobile nature of the protein. However, the simulation trajectory analysis demonstrated that the drug complex had a similar Rg trend, and no deviations were observed, denoting its stability. Hydrogen bond formation plays a vital role in providing the rigidity of a biomolecular complex. Figure 12D indicates that upon the drug binding with the M^pro complex, the hydrogen bond number did not change with respect to the control. In addition, both complexes showed similar hydrogen bond profiles over the course of the entire simulation (Figure 12D).

The binding-free energies of both systems were calculated using the MM-PBSA method. The drug complex exhibited more positive energy in the MM-PBSA calculation, indicating a more favourable binding than the control complex (Figure 13A). The control complex had a negative profile during the entire simulation period and did not exceed the drug complex in terms of positive energy. The RMSF of both complexes was also calculated to understand the degree of flexibility over the amino acid residues (Figure 13B). Most of the residues from both complexes had lower RMSF profiles, which might be responsible for the more inflexible nature of the complexes. While some deviations were observed, the RMSF value of these amino acid residues remained below 2.5 Å.

The last snapshot from the molecular dynamics trajectory was superimposed over the docked complex. The drug complex and control exhibited RMSD values of approximately 1.33 Å and 1.43 Å, respectively (Figure 14).

In silico inhibition constant (IC₅₀) calculation

The current study determined the inhibition constant of 49 towards the modelled SARS-CoV-2 M^pro. The theoretical half-maximal inhibitory concentration (IC₅₀) was calculated using AutoDock 4.0 tools, and the results of the analysis are shown in Table 7.

DISCUSSION

Because they catalyse many biochemical reactions, enzymes are crucial for numerous processes in biological systems, including the immune system. The potential role of the protease enzyme in antimicrobial responses has been previously reported in many species, such as Drosophila sp. (66). In addition, the receptor protein CD38 plays numerous immunomodulatory roles, such as those in inflammatory responses and leukocyte function (67). As mentioned above, SARS-CoV-2 M^pro plays a role in immune response interference. This study reveals the properties of Se-containing heterocyclic compounds using a series of bioinformatics tools, which could lead to the use of such compounds to increase host immune responses during SARS-CoV-2 infection.

The homology modelling technique is often used to perform the comparative modelling of protein structures because the crystal structures found in the PDB database are not always complete. In addition, compared with DNA sequences, protein sequences are more conserved. This study included homology modelling of the SARS-CoV-2 M^pro enzyme, which is essential for the viral replication cycle and polypeptide synthesis after the virus enters the host. The model generated using homology modelling most likely resembles the complete M^pro structure, with significant structural similarity (68). Furthermore, the quality of the predicted protein was verified using several standard tools, including Ramachandran statistics, Z-score and Verify 3D. In addition, we compared the findings of the Ramachandran analysis and Z-score analyses of the modelled SARS-CoV-2 M^pro with four high-resolution crystal structures from the PDB database: 6 LU7, 6 W63, 6WQF and 6XCH. Theoretically, the favoured region of a Ramachandran analysis should be a minimum of 96%; our modelled SARS-CoV-2 M^pro possessed the highest proportion of residues in the favoured region (95.1%), whereas the proportions of favoured amino acid residues for the crystal structures 6 LU7, 6 W63, 6WQF and 6XCH were 90.6%, 91.7%, 87.9% and 88.3%, respectively. However, the crystal structure 6 LU7 exhibited a higher Z-score of −7.34 than the SARS-CoV-2 M^pro model (−7.32), whereas the Z-scores for 6 W63, 6WQF and 6XCH were − 7.15, −7.11 and − 7.13, respectively. The NMA study determined that the modelled SARS-CoV-2 M^pro has a stable interaction movement. The small grooves or pockets that are responsible for the functional identification of a protein are termed the ‘active site’. In general, the active site occupies only 10–20% of the total protein but catalyses the entire enzymatic reaction (69). Therefore, the identification of the protein active site is regarded as a crucial factor for manipulating its functional properties. The dimerisation of SARS-CoV-2 M^pro has been reported to involve both Arg4 and Glu290 (70), and a previous study reported that an alanine residue at position 285 increases the catalytic activity of the protease (7). Recently, several studies have been conducted to characterise the interactions of inhibitor molecules with SARS-CoV-2 M^pro. Jin et al. showed that the inhibitor molecule N3 might serve as a potential blocker of the SARS-CoV-2 M^pro, active site, and Hung et al. reported recommended GC376 for the treatment of SARS-CoV-2 infection, demonstrating that GC376 binds with the active site of SARS-CoV-2 M^pro (43,71). In the present study, we have utilised GC376 as a positive control.

Currently, no specific therapies are available for the treatment of COVID-19 infections caused by SARS-CoV-2. However, various researchers and scientists have endeavoured to identify drug therapy options for this disease, and many studies utilising the SARS-CoV-2 M^pro have been conducted for this purpose. Several antiviral agents, including lopinavir/ritonavir, nelfinavir, remdesivir, ribavirin and favipiravir, and antimalarial drugs, such as chloroquine and hydroxychloroquine, have been used to treat COVID-19. Chloroquine and hydroxychloroquine have antiviral properties against human immunodeficiency virus (HIV), and a previous study revealed that these agents bind strongly with the crystal structure of SARS-CoV-2 M^pro (72). Although both chloroquine and hydroxychloroquine possess in vivo antiviral activity, no acute viral infection has been successfully treated by these agents in humans (73), and neither has been shown to exhibit inhibitory effect against SARS-CoV (74). Due to frequent side effects, including nausea, digestive disorders, vision impairment and heart failure, neither chloroquine nor hydroxychloroquine is currently recommended to treat SARS-CoV-2 infection (75).

The combination of lopinavir and ritonavir is approved for the treatment for AIDS, and both act as HIV protease inhibitors; ritonavir also inhibits cytochrome P450 and glycoprotein (76). A study by Muralidharan et al. revealed that a combination of lopinavir, oseltamivir and ritonavir might be effective against SARS-CoV-2 M^pro (77), and various clinical trials have examined lopinavir and ritonavir in conjunction with other drugs for the treatment of COVID-19, including arbidol, carrimycin and interferon inhalation (78). However, one study reported that the lopinavir–ritonavir combination was not associated with clinical improvements compared with standard care (79), and no positive results have been observed following oseltamivir treatment (79,80). Although another study reported that lopinavir–ritonavir combined with interferon beta-1b and ribavirin was safer than lopinavir–ritonavir alone, the risk of side effects was increased, and this finding has not yet been verified (81). The nucleoside analogue ribavirin, in combination with lopinavir–ritonavir, was found to reduce the risks of acute respiratory distress syndrome (ARDS) and death (82). However, a recent study disclosed that ribavirin required high effective concentrations for SARS-CoV-2 treatment (83). Favipiravir, an RNA-dependent RNA polymerase (RdRp) inhibitor, was shown to have more potent antiviral effects than lopinavir–ritonavir, with fewer side effects (84), but issues remain regarding favipiravir for the treatment of COVID-19 (85).

Remdesivir acts by inhibiting viral RNA polymerases (86), and previous experiments in animals showed that remdesivir could significantly decrease the viral loads in MERS-CoV–infected mouse lung tissue, improving lung function and attenuating the pathological injury to lung tissue (87). Other studies have indicated that remdesivir effectively curbs SARS-CoV-2 infections at low concentrations, and the intravenous administration of remdesivir has shown promising outcomes for the treatment of COVID-19 (14,88). However, several factors can contribute to differences in these outcomes, including the type of supportive care provided, such as medication combinations or variances in ventilation practices (89).

Se-based compounds have attracted interest in the development of novel pharmaceutical agents, including orally active Se-based antihypertensive, anticancer, antiviral, immunosuppressive and antimicrobial agents, and organoselenium compounds, which reduce oxidative tissue damage and oedema (90). In line with this, we have selected 49 potential Se-containing heterocyclic compounds for this study. In the present molecular docking simulation study, we docked these compounds with the active site of SARS-CoV-2 M^pro. A previous study reported the importance of grid-based docking (91). Our findings demonstrated that only 49, with a docking score of −6.7 kcal/mol, showed a higher affinity for SARS-CoV-2 M^pro than GC376. Intriguingly, all of the tested compounds exhibited binding affinities for residues in the SARS-CoV-2 M^pro active site in our docking analysis, and most of the compounds interacted with Cys145 and His41, forming hydrogen bond interactions with Cys145. However, only 9 formed a hydrogen bond with His41. Most of the selected compounds interacted hydrophobically with both Cys145 and His41. In comparison, GC376 only formed a hydrogen bond interaction with Cys145. Additionally, 49, 1, 3, 6, 10, 15, 26, 29, 46 and 47 formed attractive charges with Cys145. Previous studies have shown that both Cys145 and His41 form the catalytic dyad of the SARS-CoV-2 M^pro (92). In addition to interactions with Cys145 and His41, 15 compounds formed hydrogen bond interactions with Gly143; 4 compounds formed a hydrogen bond with Ser144; 39 and 43 formed a hydrogen bond with His164; 11, 13, 32, 33, 42 and GC376 formed a hydrogen bond with Glu166 residue; and 45 formed a hydrogen bond with Gln189 residue.

The substrate-binding domain of SARS-CoV-2 M^pro consists of Gly143, Ser144, His163, His164, Met165, Glu166, Leu167, Asp187, Arg188, Gln189, Thr190, Ala191 and Gln192 (93). We cross-docked the compounds that exhibited the highest docking scores with the available crystal structures of SARS-CoV-2 M^pro from the PDB database. For the current study, we selected the crystal structure of SARS-CoV-2 M^pro in complex with inhibitor N3 (PDB ID: 6 LU7) for the cross-docking analysis. The parameters were those used for the modelled SARS-CoV-2 M^pro. The cross-docking analysis using the crystal structure of SARS-CoV-2 M^pro revealed that all targeted compounds interacted with SARS-CoV-2 M^pro with high binding affinities. The positive control GC376 showed the highest docking score, followed by 49, among all targeted compounds. Intriguingly, 49 formed a hydrogen bond with Cys145 when docked with the crystal structure of SARS-CoV-2 M^pro, in contrast with the attractive charges observed for this residue when docked with the modelled SARS-CoV-2 M^pro.

The SARS-CoV-2 genome consists of 29,903 nucleotides. The 3′ terminal region primarily encodes the structural proteins, such as the N, S, membrane (M) and envelope (E) proteins, whereas the 5′ terminal encodes the two viral replicase proteins, pp1a and pp1b (94). After proteolytic cleavage, pp1a and pp1b yield 16 non-structural proteins (nsp1 to nsp16), including PL^pro and M^pro. The virus encodes other proteins, including ORF10 and nsp7b (49). These proteins can also form complexes, such as the nsp7/nsp8/Pol hetero-oligomeric complex, the spike glycoprotein/human angiotensin I-converting enzyme 2 (hACE2) hetero-oligomeric complex, the nsp7/nsp8 hetero-oligomeric complex, the nsp10/nsp14 hetero-oligomeric complex and the nsp10/nsp16 hetero-oligomeric complex (49). The SARS-CoV-2 proteome could be applied to IgG and IgM response profiling, as these antibodies have been found to bind specifically to the N and S proteins (95). Research has shown that the binding affinities between ACE2 and the SARS-CoV-2 and SARS-CoV RBDs are similar, and the RBD of SARS-CoV-2 represents a crucial region for receptor binding and could have a high potential as a target for therapeutic agents (96). Another study revealed that the pharmacological blockade of SARS-CoV-2 PL^pro not only suppressed the viral infection but also increased antiviral immunity (97). The non-structural protein nsp9 has been speculated to mediate viral replication, viral genomic reproduction and overall virulence (98). During viral infection, nsp13, nsp14 and nsp15 are recruited to double-membrane vesicles, functioning as potent interferon antagonists (99). Nsp13 also possesses NTPase and RNA helicase activity, whereas nsp15 shows RNA endonuclease activity (100,101). Therefore, we have conducted molecular docking simulations for these receptors with 49. Because the antiviral remdesivir has been shown to inhibit the RdRp protein, we also assessed the interaction between 49 and RdRp. Our study revealed that 49 interacts with the N protein, the RBD of the S protein, PL^pro, RdRp, nsp9, nsp13 and nsp15 of SARS-CoV-2. Importantly, 49 interacted with the N protein with a greater binding affinity than curcumin and mitocurcumin, which were reported in a previous study by Pal et al. (102). In addition, a previous study reported that the core structure of the RBD of the S protein, which is formed largely by three pairs, helps to stabilise the β-structure and the amino acid residues Cys336-Cys361, Cys379-Cys432 and Cys391-Cys525 from the core (96).

The current study demonstrated that 49 forms hydrogen bonds with the Trp423 and Tyr356 residues, which are located within the core and may be important for the stabilisation of the protein β-structure. We speculate that 49 might interfere with the core of the RBD of the S protein; however, co-enzyme A, the positive control used for the docking analysis of the RBD, formed more hydrogen bonds with the core amino acid residues. Despite this, 49 showed a greater binding affinity (−7.8 kcal/mol) with the RBD of SARS-CoV-2 than coenzyme-A, which showed a docking score of −6.4 kcal/mol. Compound 49 also interacted with other non-structural proteins, including PL^pro and RdRp. Specifically, our targeted compound formed a hydrogen bond with Asn109 in SARS-CoV-2 PL^pro, which has been reported to be located in the active pocket (103), in addition to hydrophobic interactions with Leu162, Glu161 and Val159, which are found in the catalytic region of PL^pro (103). Compared to 49, remdesivir showed moderately greater interaction with the RdRp of SARS-CoV-2, forming more hydrogen bonds and hydrophobic interactions. Similar to remdesivir, 49 formed a hydrogen bond with Met755, which resides in a particular motif of SARS-CoV-2 RdRp (104), and Phe753, Ser754 and Ala762 within the same motif formed hydrophobic interactions with 49. When examining nsp9, nsp1, and nsp15, 49 interacted with greater binding affinities than their respective positive controls. When interacting with the RNA replicase nsp9, 49 formed a hydrogen bond with Val41 and hydrophobic interactions with Arg39, Phe56, Ile65, Thr67 and Ile91, which have been reported to form the active pocket nsp9 (105). When interacting with the nsp13 helicase, both 49 and the positive control cangrelor possessed similar binding affinities; however, 49 showed a higher docking score than cangrelor. In addition, 49 formed hydrogen bonds with Ser289 and Arg443, and an attractive charge was visible for Lys288.

The active site residues of SARS-CoV-2 nsp13 include Lys288 and Ser289 (106), and a computational analysis revealed that several drugs, including flavin adenine dinucleotide, desmopressin, glecaprevir and rifabutin, all interacted with Arg443 (107). The binding affinity of 49 with nsp13 (−6.9 kcal/mol) was significantly greater than that of citric acid (−5.5 kcal/mol). Previously, six key amino acid residues have been reported to form the catalytic site of the uridylate-specific endoribonuclease, including His235, His250, Lys290, Thr341, Tyr343 and Ser294, and these residues are conserved in the SARS-CoV and MERS-CoV nsp15 (108). Importantly, 49 interacted with His235 residue through both a hydrogen bond and a hydrophobic interaction. In addition, 49 formed a hydrophobic interaction with Tyr343 and an attractive charge with His235.

We also evaluated the drug candidacy of the 10 compounds with the highest docking scores, which all followed Lipinski’s rule of five and Veber’s rule, with no violations. Lipinski’s rule of five suggests that compounds that violate any of the rules may struggle with permeability, absorption, of bioavailability. Compounds with lower molecular weights, fewer hydrogen bond donor or acceptor sites, and lower lipophilicity (logP) values have been associated with better permeability, faster absorption and higher bioavailability. In addition, in adherence with the parameters of Veber’s rule, the number of rotatable bonds (which predicts molecular flexibility) and the topological polar surface area predict the passive transport of a drug molecule through membranes. We also predicted the toxicological properties of the selected compounds, which showed that all compounds except 13, 26 and 30 showed compatibility with all of the toxicological parameters. In addition, the biological activities of each molecule were calculated and were compared with those of GC376. The findings revealed that all of the selected compounds showed greater ion channel inhibition and nuclear receptor ligand activities than remdesivir. The kinase inhibitor properties of 26, 44 and 49 were greater than those of GC376. Conversely, GC376 possessed stronger protease activity.

Molecular dynamics simulation studies are essential for determining the stability of a receptor–ligand complex (109). In our study, the stable nature of the docked complex was evident over a short period of time. The average RMSD value was in line with those of previous studies (110, 111), in which several druglike small molecules achieved stability in a similar time frame (111). This stability originates from the different residues and interactions present in the side chain of the complex (112). The fluctuation of the side chain was low, and only a small fluctuating region was identified, which could be due to the lack of contacts or flexible loops. Rigid structural conformations were revealed by measurements of the Rg during the simulation period (111). Superimpositions performed during the molecular dynamics simulation revealed that the post-molecular dynamics structures remained in the same cavity as the docked structure.

SARS-CoV-2 infection is characterised by not only hyperactive T cell activity but also the extreme production of inflammatory mediators, resulting in a cytokine storm (113). Traditional medicinal plants have anti-inflammatory activities capable of decreasing the level of inflammatory mediators through several mechanisms. Se has been reported to possess anti-inflammatory properties, mediated by 15-deoxy-Δ^12,14-prostaglandin J₂ (114). In addition, increased Se levels have been associated with reduced risks of malignant tumours and viral infections (23,115). Although Se can be toxic to humans, a recent study reported that high doses of Se compounds, such as sodium selenite, were well tolerated in healthy human subjects and patients (116). Based on the mechanism suggested for the Ebola virus, we speculate that sodium selenite might be used as a therapeutic agent for the prevention of coronavirus infection (117). Notably, the findings from our work agree with the findings of previous research on the selected selenocompounds, the biological activity of which has already been demonstrated. For example, ebselen (1) was shown to covalently bind the active site residue Cys145 of SARS-CoV-2 M^pro in cells, modifying the residue in an inhibitory manner (71). A previous study reported that 1 possessed the strongest inhibitory effect against SARS-CoV-2 M^pro, with an IC₅₀ of 0.67 μM (71), and analogues of 1 have been used as antiviral and antimicrobial agents (118). Research by Hsu and colleagues documented GPx-like activities for 2–5, which were moderately higher than those of their mother compound (1) (119). Viruses, such as the herpes simplex virus, have been shown to utilise the lipoxygenase (LOX) pathway to cause vascular endothelial dysfunction (120). Interestingly, 13 has been reported to inhibit LOX activity (121). Furthermore, 30 exhibited antioxidant attributes through superoxide anion scavenging, and 49 (ethaselen) is currently used as an antitumor drug in China (122). Determining the IC₅₀ represents an important step when selecting a potent inhibitor for a specific receptor. Thus, the present study also included the determination of the in silico inhibition constant of 49 against SARS-CoV-2 M^pro. The results showed that 49 exhibited an inhibition constant of 4.51 μM, which was remarkable, as the inhibition constant is directly correlated with the binding energy (64). The IC₅₀ value of GC376 has been reported as 26.4 ± 1.1 nM against SARS-CoV-2 M^pro (43). However, another study demonstrated that 49 possesses significant IC₅₀ values and inhibits the growth of several cancer cells (123). Therefore, further biological investigations may improve our understanding of the inhibitory properties of 49 against SARS-CoV-2 M^pro.

Several studies have suggested that selenocompounds may represent potent antiviral agents (124,125). Notably, selenium can be used as an adjuvant for the treatment of viral infections (126). A recent study determined that a significant concentration of selenium was found in Centella asiatica (127), which has previously been reported to be effective against viruses (128), in addition to wound healing (129) and memory-enhancing (130) and having cardioprotective (131,132), hepatoprotective (133,134), anti-filarial (135), sedative (136), anxiolytic (137), anti-tumour (138), anti-mutagenic (139), and antioxidant properties (140). Clinical trials performed to study the various activities of C. asiatica have shown efficacy against mental retardation, anxiolytic activity, venous hypertension, striae gravidarum, hepatic disorders, wound healing and vascular diseases (141).

CONCLUSIONS

Our investigations showed that Se-containing heterocyclic compounds exhibited binding affinity towards the M^pro enzyme of SARS-CoV-2, which may be associated with alterations in the enzyme’s activity. In addition, most of the compounds were found to be biologically active and possessed druglike properties. These compounds may eventually lead to the identification of potential pharmacologically active compounds, but additional in vitro and in vivo experimental validation remains necessary to confirm their efficacy against SARS-CoV-2.

Funding

The authors receive no funding from an external source.

Ahmed Rakib has completed his undergraduate and graduate studies from the Department of Pharmacy, University of Chittagong, and is now working as a research assistant in the aforementioned organization. His research interest encompasses ethnopharmacology, bioinformatics, computer-aided drug design, ion-channel pharmacology and immunoinformatics. The present address of the author is Department of Pharmaceutical Sciences, College of Pharmacy, 881 Madison Ave, Memphis, TN 38163.

Zulkar Nain is a research assistant at the Department of Biotechnology and Genetic Engineering, Islamic University, Bangladesh. His research interests are bioinformatics, functional genomics, molecular modelling as well as computer-aided drug and vaccine design.

Saad Ahmed Sami is currently completing his Bachelor of Pharmacy professional degree focused in Pharmacy from University of Chittagong. He is skilled in manuscript writing, in vivo pharmacological testing, pharmacological screening, ADME/T analysis and molecular docking simulation. His research interest mainly includes on phytopharmaceutical research and its related fields, bioinformatics and focusing on higher study and Industrial work in near future.

Shafi Mahmud is an MS student from the Department of Genetic Engineering and Biotechnology, University of Rajshahi, Bangladesh. He has interest in molecular dynamics, virtual screening, genomics and molecular biology.

Md. Ashiqul Islam is presently working as a lecturer at the Department of Pharmacy, Mawlana Bhashani Science and Technology University, Bangladesh. He completed his B. Pharm and M. Pharm from the Department of Pharmacy, University of Chittagong. His research interest lies in pharmacology, bioinformatics and computational biology.

Shahriar Ahmed has completed his B. Pharm and M. Pharm from the Department of Pharmacy, University of Chittagong. His research interests include in vivo and in vitro pharmacology, ethnopharmacology and bioinformatics.

Adnan Bin Faisul Siddiqui has successfully completed his Bachelor’s and Master’s program from the Department of Pharmacy, State University of Bangladesh. The main areas of Adnan’s research works include natural product chemistry, ethnopharmacology and molecular docking analysis.

S. M. Omar Faruque Babu currently works at the Department of Fisheries, Noakhali Science and Technology University. Omar does research in Environmental Science, Aquaculture and Agricultural Plant Science. His current project focuses on the source of metal contamination in sediment, its ecological risk and phytoremediation ability of the studied mangrove plants in ship breaking area, Bangladesh.

Payar Hossain is currently completing his Bachelor of Pharmacy professional degree focused in Pharmacy from University of Chittagong. He has completed his industrial training from Incepta Pharmaceuticals Ltd.—a leading pharmaceutical brand in Bangladesh. He is skilled in extraction methods, in vivo and in vitro pharmacological testing, and animal handling.

Asif Shahriar is a research assistant at the Department of Immunology and Microbiology, School of Medicine, University of Texas Rio Grande Valley, USA. His research interests are molecular biology, genetic engineering, virology, immunology, cancer genomics & immunoprevention and vaccinology.

Dr Firzan Nainu is currently working as an associate professor at the Faculty of Pharmacy Universitas Hasanuddin, Indonesia. Dr Nainu obtained his BSc from the Department of Pharmacy, Universitas Hasanuddin, and completed his master’s degree majoring in Biomedical Sciences from James Cook University, Australia. In 2016, he received his PhD from Graduate School of Medical Sciences, Kanazawa University, Japan. He is involved in research primarily focusing on pharmacology and immunology of infectious diseases.

Dr Talha Bin Emran has been working as an assistant professor at the Department of Pharmacy in BGC Trust University Bangladesh. Dr Emran has had his BSc honours and MS in Biochemistry from the University of Chittagong and PhD from Graduate School of Medical Sciences, Kanazawa University, Japan. His research work primarily focuses on phytochemistry, medicinal chemistry, immunology, oncology and bioinformatics.

Prof Jesus Simal-Gandara is a full professor in Nutrition and Food Science at the University of Vigo (Spain) since 1999 and vice chancellor for Internationalization in 2018. He performed research stays at the Université de Paris-Sud (Paris, France), University of Delaware (Newark, USA), Fraunhofer-Institut für Lebensmitteltechnologie und Verpackung (Munich, Germany), Central Science Laboratory (Norwich and York, UK), TNO-Voeding (Zeist, the Netherlands), Packaging Industries Research Association (Leatherhead, UK) and the Swedish Institute for Food and Biotechnology (Gothenburg, Sweden). He was 2018 and 2020 Highly Cited Researcher and has 14,850 citations in 520 papers (= 29 per paper; h-index = 65): https://www.researchgate.net/profile/Jesus_Simal-Gandara. His focus today is on the study of persistent organic pollutants from the point of view of public health (epidemiology, toxicity of mixtures, metabolites, etc.) and on the study of secondary metabolites in plant foods, exploring the molecular mechanisms that explain their activity.

References