Introduction

Hepatocellular carcinoma (HCC) is a type of primary liver cancer that develops mostly from hepatocytes worldwide (> 500,000 deaths/year)1. It is one of the most prevalent and lethal cancers in the world, having a significant impact on public health. HCC usually occurs in the setting of chronic liver illness, notably in people with cirrhosis or chronic hepatitis B or C infections, which provide an environment favorable for tumor growth2. Alcoholism, obesity, and exposure to aflatoxins are all risk factors for HCC. Early symptoms are frequently nonexistent or non-specific, making early diagnosis difficult. Patients may have symptoms such as stomach discomfort, weight loss, jaundice, and fluid collection in the belly as advanced conditions. Curative treatments for early-stage tumors may include surgical excision or liver transplantation, however advanced stages require radiofrequency ablation, transarterial chemoembolization, or systemic treatments such as targeted medicines and immunotherapies3.

The four key proteins BRAF_HUMAN, VGFR3_HUMAN, EGFR_HUMAN and UFO_HUMAN are critical in hepatocellular cancer4. BRAF is a protein kinase that is involved in the transmission of extracellular signals to the cell nucleus, consequently regulating many cellular processes. It has been estimated BRAF mutations characterize around 70% of melanoma tumors5. Vascular Endothelial Growth Factor Receptor 3, or VEGFR3_HUMAN, is a crucial protein in the processes of angiogenesis and lymphangiogenesis, its increased expression can promote tumor development and metastasis by inducing new blood and lymphatic vessels, which in turn helps provide nutrients and spread tumor cells6. Furthermore, EGFR levels are frequently observed in hepatocellular carcinoma (HCC) patients, and they have been associated with the promotion of tumor cell survival, proliferation, and motility7. Tyrosine-protein kinase receptor UFO, often known as AXL, plays a significant part in hepatocellular carcinoma (HCC). The progression of HCC has been associated with increased expression and activation of AXL, which results in a propensity for metastasis, and resistance to treatment. Nevertheless, further study is required to develop efficient medications that target this protein, particularly in the treatment of HCC8.

Resveratrol, a natural polyphenolic chemical present in many plants, especially red grapes and blueberries, has multiple features that make it a viable therapy for hepatocellular carcinoma (HCC). Its potent antioxidant and anti-inflammatory characteristics may aid in the protection of liver cells from oxidative stress and inflammation, both of which are frequently linked with HCC. Furthermore, resveratrol has antiproliferative effects in the liver, reducing the development and multiplication of cancer cells while simultaneously increasing apoptosis, which contributes to cancer cell death. The immunomodulatory properties of resveratrol may improve the immune system’s response to cancer cells, assisting in their eradication. Additionally, it can modulate several metabolic pathways and epigenetic mechanisms, which may enable cancer cells to reactivate their regular metabolic processes. Preclinical research is promising, however, there are still few human clinical studies therefore more research is needed to determine the safety and effectiveness of resveratrol as a therapeutic for HCC9.

Hepatocellular carcinoma (HCC) is a primary cause of cancer-related deaths globally. Adjuvant therapies have not been very successful in improving progression-free and overall survival, and HCC has a high probability of recurrence after treatment. Adjuvant sorafenib did not provide any appreciable benefit in the STORM study, and its use was further restricted due to the observed toxicity. There are currently trials investigating the potential of adjuvant systemic therapies, particularly immune checkpoint inhibitors (ICIs), in early-stage hepatocellular carcinoma (HCC)10,11,12,13,14,15.

The main objective of this study is to discover a potential therapeutic candidate for hepatocellular carcinoma by inhibiting proteins using a de novo drug candidate derived from resveratrol. In this study, resveratrol was fragment optimized utilizing the deep learning web server DeepFrag to enhance the phytochemical properties of resveratrol such as anti-cancerous anti-inflammatory and the results were analyzed with four targeted proteins involved in prognosis of hepatocellular carcinoma. Significant future approaches include the use of machine learning, the integration of multi-omics data, and a focus on personalized treatment. The area of study is dynamic and may see breakthroughs accelerated by collaborative research and open science initiatives, but validation through experimental testing and regulatory framework adaption will be crucial.

Material and methodology

Protein retrieval

The target proteins (ID P15056, P35916, P00533, and P30530) were retrieved from Uniprot (Universal Protein Source, (https://www.uniprot.org). The 3D structures of the selected proteins were obtained using an artificial intelligence (AI) tool called Alphafold (https://alphafold.ebi.ac.uk/), By utilizing Discovery Studio, a 3D visualisation of the protein structures was generated16.

Ligand retrieval

The literature study was performed to use phytochemical substances that contain anticancer properties and demonstrate high therapeutic potential for liver ailments. PubChem https://pubchem.ncbi.nlm.nih.gov/#query=Resveratrol) was used to obtain the structure of the Resveratrol chemical compound.

Target FDA drug

The molecular modelling and drug design lab developed the web server LigAdvisor (https://ligadvisor.unimore.it/). LigAdvisor is a flexible and intuitive tool that supports researchers in the de novo discovery and design of drugs17. The foundation of this web platform is the joint use of many ligand-based similarity techniques applied to a massive collection of publicly available chemical, biological, and structural data.

Fragment optimization

Deep Frag (http://durrantlab.com/deepfragmodel) is a deep convolutional neural network that leads ligand optimization by expanding them with chemical fragments that are highly complementary to receptors. This approach predicted how the ligand's chemical modification could enhance its binding affinity through fragment optimization by adding alternative chemical chains at a specific region of the ligand18,19.

ADME and toxicity comparison

The toxCSM online server (https://biosig.lab.uq.edu.au/toxcsm/) allows researchers to simply and accurately examine small molecule toxicity profiles, making it a useful resource. toxCSM is an efficient computational tool for analyzing and improving small molecule toxicity profiles. It outperforms previous approaches by using powerful graph-based signatures and supervised learning to generate extremely accurate prediction models20. This program provided 36 prediction models for different toxicity features, assisting in the creation of safer medications, herbicides, and pesticides for the target drug candidate designed.

Docking with pyrx

The molecular docking program with the greatest potential is Autodock Vina in PyRx (offline tool can be downloaded using the link (https://sourceforge.net/projects/pyrx/). This platform was used to determine the most suitable interacting molecules for experimental study. The molecular docking approach may be used to evaluate how tiny molecules behave at target protein binding sites and to predict the atomic interaction between small molecules and proteins21. The molecule with the lowest energy was selected for further toxicity-removal optimizations.

MD simulation

A significant computational method for examining the characteristics and behavioral patterns of molecules and materials at the atomic level is molecular dynamics (MD) simulation. The MD simulations of biological systems, such as proteins, nucleic acids, and membranes, were performed using the Desmond software package.

The MD simulation studies were carried out by utilizing the Desmond module of Schrodinger. Using its docked poses, the protein–ligand complex's dynamic behavior and stability were examined. The Protein Preparation Wizard of Maestro was utilized to preprocess the protein-ligand complex, encompassing complex optimization and minimization. The System Builder tool was used to prepare each system. The complexes were solvated using an orthorhombic box and the simple point-charge (SPC) water model. The system was neutralized with Na+/Cl− ions at a distance of 10 Å from the box’s edge. The addition of 0.15 M sodium chloride (NaCl) simulated physiological circumstances. By using the NPT ensemble, the protein complex's potential energy was reduced. For 100 ns, the molecular dynamics simulations were run at 300 K in temperature and 1 atm in pressure, with the OPLS4 force field being used for NPT creation. Before the simulation commenced, the models were relaxed. The periodic molecular systems Ewald methods22, the Coulomb interactions, 9.0 Å was the cutoff radius. The simple point charge model23 was used to explicitly define the water molecules. The temperature was controlled by the Nosé–Hoover chain coupling method24 and the pressure was controlled by the Martyna–Tuckerman–Klein chain coupling scheme with a coupling constant of 2.0 ps. Every 100 ps, the trajectories were saved for analysis, and the root mean square deviation (RMSD) of the protein complex over time was compared to confirm the stability of the simulation. For MD simulations, the predicted conformational changes from the starting structure over the course of the simulation were reported as root mean square fluctuation (RMSF) and root mean square deviation (RMSD).

Results

Proteins retrieval

The Epidermal growth factor receptor protein with Uniprot id P00533, Serine/threonine-protein kinase B-raf protein with Uniport id of P15056, Tyrosine-protein kinase receptor UFO with Uniprot id P30530 and Vascular endothelial growth factor receptor 3 protein with Uniprot id P35916 were assessed from Uniprot and further their 3D structures were retrieved from Alphafold. The 3D structures were visualized by Discovery Studio is shown in Figs. 1–4.

Figure 1
figure 1

Epidermal growth factor receptor.

Full size image
Figure 2
figure 2

Serine threonine protein kinase braf.

Full size image
Figure 3
figure 3

Tyrosine-protein kinase receptor UFO.

Full size image
Figure 4
figure 4

Vascular endothelial growth factor receptor 3.

Full size image

Phytochemical retrieval

The phytochemical Resveratrol with Puchem id of 445154 was attained from Pubchem and further the 3D structure was downloaded. The 3D structure was visualized by Discovery Studio is shown in (Fig. 5).

Figure 5
figure 5

3D structure of resverstrol from pubchem.

Full size image

FDA drug targeting

The proteins serine/threonine-protein kinase B-raf and vascular endothelial growth factor receptor 3 were given as an input for analyzing common FDA drug target against these two. The total common compounds were 4 where number of drug bank vs PDB pairs had a threshold of over 3 and common pathways were 2 including 2 Rap1 signaling pathway and focal adhesion pathway. Furthermore, the other two proteins including epidermal growth factor receptor and tyrosine-protein kinase receptor UFO were administered for common FDA drug targets. The number of total common compounds was 1 where number of drug bank vs PDB pairs had a threshold of over 7. There was no common pathway shared amongst these two proteins. The results of common FDA drugs for following proteins with MACCS similarity and ECFP4 similarity score of one is shown in the (Tables 1, 2) below. The results interpretation shows that Fostamatinib is the common FDA drug targeting all four proteins of hepatocellular carcinoma.

Table 1 The proteins serine/threonine-protein kinase B-raf and vascular endothelial growth factor receptor 3 common FDA drug targets prediction.
Full size table
Table 2 The proteins epidermal growth factor receptor and tyrosine-protein kinase receptor UFO common FDA drug targets prediction.
Full size table

Fragment optimized resveratrol

In DeepFrag online server the 0:C17 fragment, which is the cysteine residue of Resveratrol, was selected as the optimization point. The results exhibited that optimization of hydroxide ion at the selected position is evidently the right fit for optimization by scoring 1 as displayed in (Fig. 6). The Fig. 7 below show the optimization outcomes and the structures of Resveratrol before and after optimization.

Figure 6
figure 6

The different types of optimizations with varying chemical chains and their optimization score is shown.

Full size image
Figure 7
figure 7

The yellow color highlights the point 0:C17 which was optimized in resveratrol in DeepFrag webserver. The structure of resveratrol after the addition of hydroxide is shown as after optimization.

Full size image

ADME and toxicity analysis

The prediction results of ADME and toxicity generated by toxCSM server cover different aspects, which include organic, nuclear response, and genomics as are given in the (Tables 3, 4).

Table 3 Toxicity report of Fostamatinib.
Full size table
Table 4 Toxicity report of fragment optimized resveratrol (FOR).
Full size table

General molecular properties comparison

The molecular properties of weight, volume, and density including draggability properties of hydrogen bond donors, hydrogen bond acceptors and others are compared between Fostamatinib and fragment-optimized resveratrol in (Table 5) below.

Table 5 General molecular properties comparison between fostamatinib and fragment optimized resveratrol.
Full size table

Molecular drug-likeness rules

In order to evaluate the drug-likeness of Fostamatinib and Fragment Optimized Resveratrol (FOR) the molecules were assessed for the rules of Lipinski’s Rule of five, Ghose’s rules, Opera’s Notability rules, Pfizer’s rules, GSK rules and ADMETLab 2.0 Soft rules respectively. The violation of several rules are observed in case of Fostamatinib whereas, in case of Fragment Optimized Resveratrol only one rule of Oprea’s Notability rule is being violated with all other rules being fulfilled completely as shown in (Table 6) below.

Table 6 The draggability of fostamatinib and fragment optimized resveratrol were evaluated for different drug-likeness rules.
Full size table

Molecular drug likeliness score

For the estimation of the drug-likeness scores of Fostamatinib and Fragment Optimized Resveratrol the molecules were assessed for the parameters of QED scores, SAS score, Fsp3 score and NP score respectively. The Fragment Optimized Resveratrol scores higher for all aspects of its molecular properties against Fostamatinib proving it to be a better candidate for drug (Table 7).

Table 7 The molecular drug likeliness score of Fostamatinib and fragment optimized resveratrol were assessed for different properties.
Full size table

ADMET radar chart

The ADMET radar provided with the visual representation of Lipinki's rule's characteristics, such as polar surface area, logP, mass, range of atoms, range of OH, range of N or O, and range of rotatable bonds for Fostamatinib and Fragment Optimized Resveratrol. The results show that Fragment Optimized Resveratrol passed the Lipinski’s rule of five as present within the radar boundaries whereas, in case of Fostamatinib violation of Lipinki’s rule of 5 is apparent as the green lines surpassed the boundaries of optimal drug properties as shown in (Fig. 8).

Figure 8
figure 8

The ADMET radar with the representation of Lipinki’s rule’s characteristics showing fostamatinib and fragment optimized resveratrol results for following the Lipinski’s rule of 5.

Full size image

Docking analysis

The four main proteins of hepatocellular carcinoma were docked with Resveratrol and Fragment fragment-optimized resveratrol as shown in (Figs. 9–12). The energies were optimized after docking of all four proteins with Fragment Optimized Resveratrol proving its better efficiency as shown in (Table 8) below.

Figure 9
figure 9

Docking complex of EGFR_HUMAN with FOR.

Full size image
Figure 10
figure 10

Docking complex of VGFR3_HUMAN with FOR.

Full size image
Figure 11
figure 11

Docking complex of UFO_HUMAN with FOR.

Full size image
Figure 12
figure 12

Docking complex of BRAF_HUMAN with FOR.

Full size image
Table 8 The binding affinities of resveratrol and fragment optimized resveratrol with respective proteins.
Full size table

Molecular dynamic simulation evaluation

Molecular dynamics simulations were performed on the top hits containing high binding energies. Over the simulation period, the projected conformational changes from the initial structure were presented in terms of root mean square deviation (RMSD). Moreover, structural stability, atomic mobility, and residue flexibility at times of interaction of protein-hit were expressed with root mean square fluctuation (RMSF) values. The peaks of RMSF graph represent the fluctuation portion of the protein through the simulation. The N- and C-terminal show greater fluctuations overall. Alpha helices and beta strands show high fluctuation as well. All protein frames are first aligned on the reference frame backbone, and then the RMSD is calculated based on the atom selection. Monitoring the RMSD of the protein can give insights into its structural conformation throughout the simulation. The RMSD of FOR and VGFR3_HUMAN complex showed small deviation from almost 2 to 8.5 Å till 10 ns and then it was stable at almost 9 to 10 Å from 10 to 30 ns. After a peak at 16 Å on 40 ns the system remained equilibrated till 90 ns with one peak at 16 Å at almost 95 Å towards the 100 ns frame (Fig. 13). After looking the trajectories, it was found that the both systems were stable and ligands remained inside the binding pockets and made important interactions and the backbones were consistent. Similarly, estimated RMSF values of up to 4 Å indicated high stability of the complex whereas fluctuation changes at almost 5 Å is the highest where ligand made interaction with the receptor shown in green lines (Fig. 14).

Figure 13
figure 13

Root mean square deviation plot of disentanglement docked complex of FOR and VGFR3_HUMAN.

Full size image
Figure 14
figure 14

Root mean square fluctuation plot of disentanglement docked complex of FOR and VGFR3_HUMAN.

Full size image

To determine the ligand's effect on whole protein, the Disentangle docked complex of FOR and VGFR3_HUMAN was investigated. During the 100 ns simulation, six properties were taken into consideration at in order to demonstrate the stabilities of the chosen ligands in the binding pocket. Ligand RMSD, or root mean square deviation, is the measurement of a ligand's deviation from the reference conformation, which is usually the first frame, which is considered to be time t = 0. (2) Radius of gyration (rGyr): This quantity is equal to a ligand's primary moment of inertia and is used to gauge how "extended" a ligand is; The quantity of internal hydrogen bonds (HB) in a ligand molecule is known as the intramolecular hydrogen bond (intraHB). (4) Molecular surface area (MolSA): A probe radius of 1.4 Å was used to compute the molecular surface; (5) Polar surface area (PSA): Solvent accessible surface area in a molecule given only by oxygen and nitrogen atoms; (6) Solvent accessible surface area (SASA): Molecular surface area of accessible by a water molecule (Fig. 15).

Figure 15
figure 15

Variation in the ligand’s properties in disentanglement docked complex of FOR and VGFR3_HUMAN with respect to time during the course of 100 ns simulation.

Full size image

Ligand RMSD

As shown in Fig. 15 the RMSD of the ligand was 0.5 Å till 10 ns and after that it was stabilized at 0.75 Å till 100 ns.

Radius of gyration (rGyr)

The radius of gyration measures the ‘extendedness’ of a ligand, and is equivalent to its principal moment of inertia. The rGyr value of ligand was from 4.16 to 4.3 Å at the start and remained constant at almost 4.20 Å till the end. The constant values indicate the steady behaviour in (Fig. 15).

Intramolecular hydrogen bonds (intraHB)

The number of internal hydrogen bonds (HB) within a ligand molecule showed by intraHB are not found in this particular simulation of the FOR and VGFR3_HUMAN complex as shown in (Fig. 15).

Molecular surface area (MolSA)

Molecular surface is the calculation with 1.4 Å value of probe radius. This value is equivalent to a van der Waals surface area. The MolSA value was 242 to 244 Å2 throughout remining stabilized till the time of 100 ns as shown in (Fig. 15)

Solvent accessible surface area (SASA)

SASA is the Surface area of a molecule accessible by a water molecule. The value of SASA was almost 600 Å2 till the 100 ns of simulation as shown in (Fig. 15). No unusual fluctuations in the graph were observed throughout the simulation representing the steady behavior of ligand.

Polar surface area (PSA)

The PSA is Solvent accessible surface area in a molecule contributed only by oxygen and nitrogen atoms. The PSA value was 192 to 196 Å2 till 100 ns timeframe depicting complete stability throughout as shown in (Fig. 15). Several energy components inside a system can be analyzed to get insight into prevailing forces and reach simulated equilibrium. An easy-to-use interface is provided by the Trajectory Plot tool for configuring and performing these kinds of analysis. It is necessary to have access to a Linux-based host computer with a compatible GPU card in order to perform energy analysis. To investigate the ligand’s internal energy, a Linux system with two GPU cards—the RTX 3070 and 3080—was used in this work. Correlating energy fluctuations with conformational changes and examining the energy exchange between a protein and receptor over the course of a 100 ns simulation were the objectives of this investigation. With an indicated value of −21514 cal/mol, the energy breakdown covering Coulomb (electrostatic) and van der Waals (vdW) interactions demonstrated a significant binding strength (Fig. 16).

Figure 16
figure 16

Receptor-ligand complex energy analysis throughout 100 ns simulation.

Full size image

MM/GBSA binding free energy of complex

The Desmond's thermal_mmgbsa.py script and Molecular Mechanics Generalised Born Surface Area (MM/GBSA) were used to calculate the binding free energies. The following formula is used to compute the MM/GBSA free energy of binding (ΔGbind) for the docked postures and Desmond trajectories:

$$\Delta {\text{GBind }} = \, \Delta {\text{GComplex }} - \, \Delta {\text{GLigand }}{-} \, \Delta {\text{GReceptor}}$$

where the variables ΔGComplex, ΔGLigand, and ΔGReceptor in MM/GBSA reflect the energy estimations of the optimised complex (complex), optimised free ligand (ligand), and optimised free receptor (receptor).

MM/GBSA is unquestionably one of the most popular methods for estimating binding free energy because it is more accurate than the majority of molecular docking score systems. MM/GBSA is widely used in bimolecular research describing, among other things, protein folding, protein-ligand binding, and protein-protein interaction since it uses less processing power than other chemical free energy scoring techniques. One of the most crucial problems in bimolecular investigations is the precise calculation of the free energy, which drives molecular motions. The solvation free energy difference upon binding and the interaction energy between the ligand and receptor complex can be used to compute the binding free energy using the MM/GBSA method.

Throughout the post-MM-GBSA computation of 100 ns MD data of both hit ligands identified by the dynamic’s investigations, a total of 1000 frames were processed and examined. Every 10 ns, we performed binding free energy estimations using 11 recorded images. The mean ΔG was discovered to be −70.80 kcal/mol. The energies were displayed in (Fig. 17) every 10 ns.

Figure 17
figure 17

MM/GBSA approach to calculate Binding free energy of complex.

Full size image

Discussion

Globally, hepatocellular carcinoma (HCC) is the most prevalent primary liver cancer and the major cause of cancer-related mortality. HCC ranks as the ninth most common cause of cancer-related fatalities in the US. Despite advances in screening, prevention, and diagnostic and treatment technologies, the incidence and death rates continue to rise25. The four important proteins BRAF_HUMAN, VGFR3_HUMAN, EGFR_HUMAN, and UFO_HUMAN are involved in cell growth, proliferation, and differentiation in hepatocellular carcinoma. Increased expression and activation of these proteins leads to increased aggressiveness, propensity for metastasis, and treatment resistance26.

Immuno checkpoint inhibitors (ICIs) and their combinations have revolutionised the treatment of hepatocellular carcinoma (HCC). In comparison to sorafenib, atezolizumab plus bevacizumab considerably improved survival in the IMbrave150 study, which established it as the new standard. Trials combining atezolizumab and cabozantinib as well as sintilimab and a biosimilar bevacizumab are also demonstrating promise.

Alternative and complementary medicine, including plant-based medications with antioxidant and neuroprotective properties, has gained popularity in recent years. Resveratrol, a polyphenolic molecule present in many plant species, has emerged as a promising nutraceutical with potential therapeutic use in neuropsychiatric, cardiometabolic, and cancer illnesses, as well as aging. The effects of resveratrol are biphasic and dose-dependent27. It functions as an apoptotic agent for cancer prevention at high doses in the range, and provides cardioprotection at relatively low levels by acting as an anti-apoptotic agent. Resveratrol’s possible toxicity was assessed by toxicity predictions, and the findings showed that resveratrol is adverse to the kidneys. Despite this, we were able to create resveratrol fragments from the Deepfrag bioinformatics tool that are non-toxic and produce more effective outcomes28.

Fostamatinib, also known as Tavalisse, is an oral spleen tyrosine kinase (SYK) inhibitor used to treat adult patients with persistent immune thrombocytopenia (ITP) who have not reacted well to conventional medications. This chemical controls immunological responses by inhibiting SYK, a critical protein kinase implicated in a variety of immune and inflammatory processes, making it a prospective treatment option for hepatocellular carcinoma. While oral administration is convenient, Fostamatinib have no side effects. It targets all four proteins (braf_human, vgfr3_human, egfr_human, and ufo_human) of hepatocellular carcinoma and has no adverse effects in the human body29. Fostamatinib Violate no. 2 the Lipinski's RO5, no. 2 the Ghose's Rules, no. 1 the Pfizer's Rules, no. 1 the GSK Rules, and no. 4 the ADMETLab 2.0 Rules. However, in this work, we created a de novo candidate by fragment optimization of resveratrol, which produced more accurate and superior outcomes with only one violation, as shown in (Table 6). Since there are so few medications in clinical trials, this de novo-developed drug is less toxic and has superior drug likeness, and this designed drug with different analyses demonstrated it to be the best candidate30,31,32.

This study’s reliance on bioinformatics tools obliges authentication through in vitro and in vivo experiments to confirm the efficacy and safety of the resveratrol-derived compound. Moreover, generalizability is limited due to the small sample size and diversity, potential long-term side effects, and unexamined drug interactions. Future research should address the limitations to improve the applicability of the findings.

Conclusion

The most frequent kind of primary liver cancer is hepatocellular carcinoma (HCC). Hepatocellular carcinoma is more common in persons who have chronic liver disorders, such as cirrhosis caused by hepatitis B or hepatitis. Resveratrol is a natural polyphenol present in some plants that has been studied for its possible health advantages, including antioxidant and anti-inflammatory effects. However, in this study, we developed a de novo candidate by fragment optimisation of resveratrol, which generated more accurate and superior results and would deem a better pharmacological candidate than Fostamatinib. In future in-vitro and in-vivo studies can be performed.