An Evidence of Drug Repurposing for COVID-19 Pandemic Based on In silico Investigation from Phenolic Derivatives of Silybum Marianum Against SARS-Cov-2 Proteins

Swaraj Mohanty ¹, Soumya Lipsa Rath ², Poornima Sharma ¹, Yasmin Ahmad ^1,*

1. Pathophysiology and Disruptive Technology (PDT) Division, Defence Institute of Physiology & Allied Sciences (DIPAS), Defence Research & Development Organization (DRDO), Lucknow Road, Timarpur, Delhi-110054, India

2. Department of Biotechnology, National Institute of Technology (NIT), Warangal, Telengana-506004, India

* Correspondence: Yasmin Ahmad

Academic Editor: Xiangzi Han

Special Issue: Oxygen Transport Physiology and COVID at High Altitude

Received: January 11, 2023 | Accepted: June 19, 2023 | Published: July 04, 2023

OBM Genetics 2023, Volume 7, Issue 3, doi:10.21926/obm.genet.2303186

Recommended citation: Mohanty S, Lipsa Rath S, Sharma P, Ahmad Y. An Evidence of Drug Repurposing for COVID-19 Pandemic Based on In silico Investigation from Phenolic Derivatives of Silybum Marianum Against SARS-Cov-2 Proteins. OBM Genetics 2023; 7(3): 186; doi:10.21926/obm.genet.2303186.

© 2023 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

In late 2019, the outbreak of COVID-19 completely changed the global healthcare scenario. This highly communicable disease created an emergency healthcare condition across the globe. The countries were also affected during this global pandemic due to poor hospital management conditions and economic crashes. Nevertheless, considering the situation's urgency, companies worldwide worked towards generating vaccines and alternative therapeutics to eradicate the disease globally. However, this objective could not be achieved completely due to the emergence of new variants of the virus and the increasing drug resistance and escape strategies of the virus.

Developing a new drug or a vaccine is tedious, from drug discovery to in-vitro screening, preclinical and clinical trials. Gaining FDA approval requires a minimum of 12 to 15 years, which is too long considering the mortality rate of COVID-19. An alternative strategy that can boost the process is to repurpose existing drugs. Repurposing FDA-approved drugs drastically reduces the time and costs incurred in manufacturing. Major flavonolignans such as silibinin, isosilybin, silychristin and silidianin from Silybum marianum are most prominently known for their medicinal properties to cure different disease states of the liver like chronic liver cirrhosis [1], fibrosis [2], necrosis [3,4] and hepatocellular carcinoma [5,6]. This compound is also well known for its anti-inflammatory, antioxidant, anti-carcinogenic and anti-mutagenic properties. Silymarin on the rat as an animal model demonstrated that it helps in the prevention of free radical generation of reactive oxygen species (ROS), triggers the antioxidant defense system, improves the anti-inflammatory response and activates the key genes that act as compensatory adaptive vascular response during the hypoxic condition [1,7]. As silymarin is in the clinical trial pipeline and its less toxic properties is a unique characteristic that can ensure it as an effective drug molecule for a different disease, we evaluated the efficacy of the molecule towards main protease, spike receptor binding domain and RNA binding domain of nucleocapsid phosphoprotein enzymes of SARS-CoV-2 [8]. The COVID-19 disease affects the endothelial lining of the lungs and the disease progression also relates to physiological changes as in the case of high-altitude hypoxia illness [9]. Hence, we want to check the possibility of silibinin (a major constituent of silymarin) as an effective inhibitor.

The SARS-CoV-2 genome consists of 12 open reading frames (ORFs), 9 transcription regulatory sequences, 9 conserved leader sequences and 2 untranslated regions (UTRs) [10] (Figure 1). The multiple sequence alignment data of the SARS-CoV-1 & SARS-CoV-2 shows the sequence similarities, and it is believed that 10-28 nucleotides in the 5' UTR interacts with the non-structural protein1 (Nsp1). Targeting the major viral protein to stop the progression of pathogenesis needs a clear understanding of the viral genome and its translational units that interact with different host machinery, starting from the viral entry till the hijack of the cell immunity [11]. Among the four structural proteins of the SARS-CoV-2 virus we have targeted the binding sites of major proteins, i.e., the RNA binding domain of the nucleocapsid of phosphoprotein and the receptor binding domain of spike protein. From the virology aspects main proteases play a major role in viral proprotein maturation and assembly of other components to produce a whole virus after the translational event inside the host cell. We also targeted the cysteine-like protease, i.e., the main protease of the ORF1a region [8,12]. Hence investigation to find the major inhibitory molecule that can actively bind and block the enzymatic activity of the major proteins of SARS-CoV-2 is necessary to break the chain of the disease transmission from the pharmacological point of cure.

Figure 1 The genome organization of SARS-CoV-2 and its translational regions.

2. Materials and Methods

2.1 Molecular Docking and ADME/Tox Analysis

The 3D structural data files (SDFs) of ligands (CID: 31553, 3085830, 441764 & 1982272) of Silybum marianum were taken from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) and checked for the pharmacokinetics and toxicity properties with ADMElab 2.0 (https://admetmesh.scbdd.com/). This ADMElab 2.0 tool is considered to efficiently calculate and predict 17 physiochemical, 13 medicinal chemistry properties, 23 different ADME, and 27 toxicity endpoints along with the 8 toxicophore rules.

The crystallographic structures of the proteins of SARS-CoV-2 were imported from the protein data bank (PDB ID: 6LU7, 6VYO & 7BZ5). The ligands were transformed into PDBQT file format by using Open Babel (Version 3.1.1) tool [13]. The refinement of protein structures was done with AutoDock Vina (version 1.2.0) tool developed at The Scripps Research Institute [14]. The protein structures were imported, and the water molecules present in the crystal structure data were removed. Interacting small molecules were also removed from the structure, and missing residues and polar hydrogen molecules were appended. The refined structures were subjected to energy minimization and then a grid box with a dimension of X, Y & Z coordinates = 126.0 and spacing 0.375 Ang was used for blind docking. The center of the grid was fixed as: -center_x = -26.286, center_y = 12.608 & center_z = 58.965 for 6LU7, center_x = -14.253, center_y = 41.224 & center_z = 14.164 for 6VYO and center_x = -72.527, center_y = -29.929 & center_z = 11.449 for 7BZ5. Gasteiger charges were added to the protein before docking. All rotatable bonds of ligands were kept flexible while the protein remained rigid during the docking.

2.2 Molecular Dynamics Simulation

Based on the strong docking interaction profile and the lesser toxicity for the silibinin compound in the ADMET results, we focused on studying the molecular evolution and the conformational changes of macromolecules (PDB ID: 6LU7, 6VYO & 7BZ5) with the ligand molecule during the 100 ns molecular dynamics simulation. Molecular Dynamics (MD) Simulation was performed by using GROMACS (version 2020.4) MD simulation software package [15,16]. Charmm36-jul2020 force field parameters for protein were used for the study [17,18]. The topology and parameters for ligand molecules were generated using the CHARMM General Force Field Server [19]. The complexes were placed inside a box filled with TIP3P water molecules and neutralizing ions. The box edge was kept at 1 nm from the protein-ligand complex to avoid calculation artifacts. Initial minimization was performed using steepest descent and conjugate gradient algorithms followed by equilibration at an isothermal-isobaric ensemble. A time step of 2 fs was considered during the simulation and A trajectory of 100 ns was generated to obtain data for analysis. Pymol [20,21] and BIOVIA Discovery studio visualizer [22] were used to analyze structures generated during the study.

3. Results and Discussions

3.1 Description of Silibinin & Its Derivatives

3.1.1 ADME/Tox Properties

The absorption of any oral drug molecule takes place in the intestinal cells before releasing into the systematic circulation with an active and passive diffusion process. Hence different in vitro models have been developed for the permeability study and oral bioavailability. Likewise, the distribution of the drug after administration has also been calculated along with the concentration-to-volume ratio, plasma protein binding, and the side effect in the central nervous system (CNS). The metabolism of any drug molecule occurs in the liver using oxidative reaction and conjugative reaction; isozymes of the cytochrome P450 family have been studied. Clearance of the drug molecule, its half-life and toxicity are also crucial parameter in understanding the pharmacokinetic of drug [23]. The four compounds show good medicinal properties as they satisfy the Lipinski rule, Pfizer rule and Golden triangle rule with the silibinin and isosilybin showing medium apparent permeability coefficient(P_app) to estimate the effect of blood-brain barrier(BBB) by considering the Maidin-Darby Canine Kidney cells(MDCK) as an in vitro model. With a moderate P-glycoprotein(P-gp) inhibitor property and excellent P-glycoprotein substrate forming ability of Silibinin and isosilybin is known to protect the body as it maintains the removal of drugs from the kidneys and liver and maintains the integrity of BBB. The metabolism of the intake drug plays a major role in keeping the homeostasis of major organelle of an individual with complex physiological interaction; the very minimal(<1) values of CYP34A substrate and inhibitors from this drug profiling explicitly cite that these are the major molecule accounts for 30-50% of drug metabolites. Silibinin & isosilybin shows a very lower dose intake as per the mmol/kg- bw/day as per the calculation by Food and Drug Administration maximum daily dose(FDAMDD) and lower rat oral acute toxicity which supports these two isoforms to be with an excellent F_20% oral bioavailability. With a higher acceptance to the Toxicology in the 21^st century(Tox21) methodology by database screening and excellent scoring of nuclear androgen receptor(NR-AR), NR-AR binding to ligand(NR-AR-LBD), nuclear estrogen receptor(NR-ER) and NR-ER ligand binding affinity these two molecules needs further attention to the pharmacokinetics studies. The effect of the chemical molecules can be well understood from the values (Table 1) and the explanation of different properties are supplied as a supplementary datasheet (Supplementary datasheet).

Table 1 ADMET Property of Ligands.

3.1.2 Other Characterization

The physicochemical and medicinal chemistry properties of ligand molecules are properties satisfying and the molecules can be called potent drug candidates and need further clinical and non-clinical studies to fall under drug development pipelines. The drug-likeness properties of the ligand molecules can better be understood from the bioavailability radar. The parameters (like molecular weight, van der Walls volume, density, number of hydrogen bond acceptors, number of hydrogen bond donors, number of rotatable bonds, number of rings, number of atoms in the biggest ring, number of heteroatoms, formal charge, number of rigid atoms, flexibility, number of stereocenters, topological polar surface area, the logarithm of aqueous solubility value, logarithm of the n-octanol distribution co-efficient) were studied and represented (Table 2) (Figure 2).

Table 2 Physicochemical & Medicinal Property of Ligands.

Figure 2 The image represents the bioavailability radar for the ligand molecules silibinin (a), isosilybin (b), silychristin (c) and silidianin (d).

3.2 Silibininas A Potential Drug Molecule

3.2.1 Molecular Docking Results

The molecular docking of different proteins of SARS-CoV-2 with flavonoid compounds shows good binding scores and interaction with the active residue of proteins with the binding site of ligand molecules (Table 3). The data was analyzed from the docked complex with the help of PLIP [24] and BIOVIA discovery studio visualizer [22]. Silibinin molecule shows higher negative binding affinity with the three major macromolecules and the nearest bonding forming ability with major amino acids present in the pocket region of the macromolecules. The electrostatic/hydrophobic interaction of essential amino acids like methionine, tryptophan, and phenylalanine in the case of main protease, spike protein and the non-essential amino acids like proline with the nucleocapsid phosphoprotein respectively. The 3D and 2D interactions of the complexes have been represented (Figure 3, Figure 4, Figure 5).

Table 3 Molecular Docking Results.

Figure 3 Interaction of main protease (6LU7) with different ligand molecules (silibinin, isosilybin, silychristin & silidianin) in 3D diagrammatic view (a) and 2D schematic view (b).

Figure 4 Interaction of RNA binding domain of nucleocapsid phosphoprotein (6VYO) with different ligand molecules (silibinin, isosilybin, silychristin & silidianin) in 3D diagrammatic view (a) and 2D schematic. View (b).

Figure 5 Interaction of spike receptor binding domain (7BZ5) with different ligand molecules (Silibinin, isosilybin, silychristin & silidianin) in 3D diagrammatic view (a) and 2D schematic view (b).

3.2.2 Bioactivity Score Analysis

The four lead molecules silibinin, isosilybin, silychristin and silidianin were subjected to bioactivity score analysis based on the parameters like G protein-coupled receptors (GPCR), Ion channel modulator (ICM), Nuclear receptor ligand (NRL) and inhibitory enzymes (Protease and Kinase). The scores greater than 0.00 considered highly active; values ranging from -0.50 to 0.00 were moderately active and those less than -0.50 were inactive [25]. Based on the analysis score predicted silibinin and isosilybin were as effective due to their higher enzymatic activity inhibitor and hence can be taken forward into further drug development (Table 4).

Table 4 Bioactivity Score of Ligands.

3.3 Stability of Silibinin in Protein Pockets by MD Simulations

A molecular dynamic simulation of 100 ns was analysed using the different GROMACS modules with the pre-set algorithm and the results have been explained.

3.3.1 Protein Characterization

Root Mean Square Deviation (RMSD). The root means square deviation (RMSD) trajectory for a simulation run for three macromolecules with silibinin molecule shows the overall stability of the protein-ligand complex during the binding with active site amino acid residues of the protein molecules. The protein molecules are globular so the acceptable deviation ranges within a difference of 1-3 Å. From the C-α RMSD plot of the main protease (PDB ID: 6LU7) we can observe the trajectory lies within the range with very minimal deviation and after 80 ns up to 100 ns both the protein and the ligand RMSD shows a good binding affinity as in equilibrium and without any fluctuation. Likewise, if we look into the trajectory of the nucleocapsid protein (PDB ID: 6VYO) we can clearly understand that the protein molecule is very stable. The fluctuations are within 2.5 Å and most importantly during the complete 100 ns simulation, and the molecular ligand affinity shows a flattened trajectory from 78 ns up to 97 ns. And when we considered the RMSD trajectory for spike receptor binding domain (PDB ID: 7BZ5) with the ligand complex we observed the protein molecule gradually came to equilibrium in fluctuation from 68 ns up to 100 ns and surprisingly the ligand fluctuation is within 0.5 Å throughout the 100 ns simulation run (Figure 6a).

Figure 6 (a) Root mean square deviation (RMSD) trajectories of macromolecules (6LU7, 6VYO & 7BZ5) (BLACK) with silibinin molecule (RED) and (b) Root mean square fluctuation (RMSF) trajectories of macromolecules (6LU7, 6VYO & 7BZ5) during 100 ns molecular dynamics simulation.

Root Mean Square Fluctuation (RMSF). The root mean square fluctuation (RMSF) of protein and ligand molecules gives the idea of each residue present in molecular structure. This calculation generates a trajectory that denotes the flexibility of individual amino acid residues and atoms during the simulation process. When we observed the RMSF trajectory of the macromolecules we could observe the flexibility of all 306 amino acid residues of the main protease (6LU7), 128 amino acid residues of the nucleocapsid protein domain (6VYO) and 229 amino acid residues of the spike protein (7BZ5) (Figure 6b). The first 3 amino acid residues of the 6LU7 showed a high fluctuation of 6 Å. However, the fluctuation reduced to within a difference of range within 3 Å for other residues throughout the simulation run. In the case of the 6VYO protein, the fluctuation is minimal and ranges between 0.5-2.5 Å indicating no more stretching of the bonds formed between the active site residues of the complex. The trajectory of the 7BZ5 molecule also showed minimal fluctuation for a maximum number of residues during the simulation process.

Principal Component Analysis (PCA). The protein-ligand interaction results in a thousand possible poses during the simulation process and it is very difficult to analyze every pose without a statistical tool. PCA is used as a mathematical tool to detect the correlation between a large set of datasets; its biological applications are to detect the flexible regions in a protein that hinder the equilibrium state of the protein. This PCA can help integrate physical models of protein motions after removing atoms’ translational and rotational movements when interacting with drug molecules. The distribution of atoms during the simulation process has been represented in the form of dot Cartesian trajectory coordinates from most eigenvector values. From our result interpretation we can predict that the interaction of drug molecules with nucleocapsid phosphoprotein and spike proteins shows very minimal changes whereas the structure of the main protease has been distorted with a large variation. The eigenvalues and eigenvectors of the covariance matrix were diagonalized with the first two principal components, i.e., PC1 & PC2 and the first 25 eigenvectors were considered (Figure 7). We conclude the silibinin drug shows maximum effectiveness against the main protease (i.e., 6LU7) and can be used as a target against the maturation of viral accessory polyproteins inside the host cells and hence be helpful for the retardation of viral proliferation [26,27,28]. When we observed the PCA plot for nucleocapsid phosphoprotein and spike glycoprotein we can state that residues of proteins showed less movement during the simulation process and hence preferred to remain in a natural state.

Figure 7 PCA of different macromolecules (6LU7, 6VYO & 7BZ5) and eigen values of the covariance matrix during 100 ns simulation process.

Secondary Structure (SS). During the protein-ligand interaction, secondary structural changes occur in protein molecules along with the evolution of time using the DSSP command line tool of GROMACS. The changes in the elements like helix and beta-sheets were thoroughly observed for all the residues of the three proteins and ligand binding complexes (Figure 8). A very minute fluctuation occurred in the turns, α-helix and 3-10 helix of 6LU7 (main protease). In the case of 6VYO (nucleocapsid phosphoprotein) we can visualize that the turns, β-sheets, β-turns, 3-10 helix and coils show a minimal alteration in natural structures. Similarly, if we consider the 7BZ5 (spike protein) we have observed the distortion in the second segment as compared to the natural structure in the turns, α-helix, π-helix, 3-10 helix whereas β-sheets and β-turns show minimal fluctuations. Since the structural integrity of protein depends on the backbones of protein, a protein’s secondary structure plays a vital role in protein folding and misfolding [29]. In the case of the intervention of drug molecules the distortion of protein structure occurs and protein functions and bioactivity get disturbed [30,31]. The in silico interaction of the silibinin molecule with the main protease and the spike protein shows promising results in distortion of the secondary structures at different backbone residues and hence can affect the protein conformational changes leading to protein functionality [32].

Figure 8 Secondary structure of macromolecules (6LU7, 6VYO & 7BZ5) during 100 ns MD simulation process.

3.3.2 Ligand Characterization

Root Mean Square Fluctuation (RMSF). We have also analysed ligand flexibility of the molecules of silibinin compound during a simulation run the trajectory for all atoms and observed deviation falls within the difference of maximum 3 Å for all three different molecular binding sets and represented below (Figure 9). The more the atoms’, the more the atoms’ flexibility to bind with active site residues of the protein molecules leads to stronger molecular interaction among protein and ligand molecules [33]. It is being observed that all three proteins formed efficient binding interaction with silibinin molecules during the 100 ns MD simulation process and hence can be considered as a good inhibitory molecule [34].

Figure 9 RMSF of silibinin ligand with different macromolecules (6LU7, 6VYO & 7BZ5) during 100 ns MD simulation.

Radius of Gyration (Rg). The radius of gyration is defined as the distribution of atoms in a protein molecule around its axis. During the binding of a lead compound with protein the conformational changes in structure are observed as the compactness of the molecule gets disturbed due to different binding forces. The lesser deviations in values from the central axis during the simulation process more the structural integrity of the molecules is preserved [35]. Hence, the Rg plotted below shows much less fluctuation within a difference in the 0.5-1 Å (Figure 10a). This also predicts compactness of macromolecules is not getting disturbed when binding with silibinin inside the in silico environment [36].

Figure 10 Radius of Gyration (a) and solvent accessible surface area (b) of different macromolecules (6LU7, 6VYO & 7BZ5) during 100 ns MD simulation process.

Solvent Accessible Surface Area (SASA). The protein-ligand binding is a solvent-substitution process in which the protein gets unfolded to provide a surface area on the active site residues of the proteins when exposed to a solvent system. After molecular dynamics simulation, it is highly necessary to compute the SASA which shows how efficiently the computational system can mimic the intracellular physiological environment [37]. The low fluctuation in the SASA value during the simulation process suggests the very good inhibitory action of the molecule against the protein targets [38,39]. While the binding of the molecules with the ligand molecule parts of the proteins are buried inside the solvent and few portions are exposed to the solvent; the trajectory below represents that the fluctuation of the SASA values is within a different range of 150 Ų and are the surface that is exposed into the solvent system (Figure 10b).

3.3.3 Protein-Ligand Interaction Profile

The molecular dynamics simulation of different proteins of SARS-CoV-2 with silibinin compounds shows good interaction during the active pocket residue of proteins with the binding site of ligand molecules (Table 5). The last frame complex during the simulation process was analyzed with the help of PLIP [24] and BIOVIA Discovery Studio Visualizer [22]. The 3D interaction of the complexes and the 2D interactions of complexes have been represented (Figure 11). We conclude that among the interaction of three macromolecules the maximum number of hydrogen bonding with bond orders less than 3.5 Å was seen between active site residues of the main protease and ligand molecule. Whereas several hydrophobic interactions have been found between ligand and spike protein, since the complete process uses water as a solvent system, we can elucidate that the hydrogen bonding is the strongest interaction and more efficient attachment of drug candidate to the target.

Table 5 PLIP Molecular Dynamics Simulation Results.

Figure 11 Protein-ligand interaction profile (PLIP) of silibinin with different macromolecules (6LU7, 6VYO & 7BZ5) in 3D (a) & 2D (b) schematic view representing different bonds and bond length after 100 ns MD simulation last frame.

4. Conclusions

It would be a great aid to repurpose this medical emergency the surge in cases due to COVID-19 worldwide researchers are striving to find an efficient drug candidate to alleviate the disease progression and treatments. These study outcomes can facilitate the silibinin molecule as a promising molecule for binding to the main protease and inhibit the formation of viral accessory assembly and replication of the SARS-CoV-2. And electrostatic interaction of the silibinin molecule with the spike protein of the virus can impede the interaction of the ACE2 receptor of the host by blocking the receptor binding domain (RBD) of the spike protein. Since COVID-19 shows critical pulmonary endothelial dysfunction by clogging atrial and venous fluid flow, the reported anti-inflammatory and anti-thrombotic properties of silibinin molecule in treating hepatic diseases will be relaxing for complications during the disease progression. The reported study on the size reduction and therapeutics intervention in animal models and its cell permeability data and the pharmacovigilance study are available as this drug has been used for many years in patients suffering from liver diseases. The candidate drug silibinin should be further tested in the COVID-19 drug repurposing pipeline and if necessary, the therapy may include other existing therapeutics based on the patient’s clinical condition. We believe this silibinin is a lifesaving drug in case of COVID-19 management as a multitarget drug, without any harmful side effects and investment of huge financial cost, time and effort.

Author Contributions

YA conceptualized the study. SM performed the ADMET, molecular docking (MD) and wrote the entire manuscript with the necessary tables and sketched figures for convenient representation. SLR set the environment and carried out the molecular dynamics simulation (MDS) by using the computational resource facility of National Institute of Technology (NIT)-Warangal. SM and SLR equally contributed in this manuscript hence share the equal authorship. SM and PS checked the language and formatting of the manuscript. YA and SM critically evaluated the manuscript.

Funding

SM is a recipient of DRDO-SRF fellowship and PS is a recipient of DRDO-JRF fellowship. The authors have not received any special funding for this work apart from this institutional financial assistance.

Competing Interests

Authors state no conflict of interests.