Geometric Deep learning Prioritization and Validation of Cannabis Phytochemicals as Anti-HCV Non-nucleoside Direct-acting Inhibitors.

Abstract

Introduction:

The rate of acute hepatitis C increased by 7% between 2020 and 2021, after the number of cases doubled between 2014 and 2020. With the current adoption of pan-genotypic HCV therapy, there is a need for improved availability and accessibility of this therapy. However, double and triple DAA-resistant variants have been identified in genotypes 1 and 5 with resistance-associated amino acid substitutions (RAASs) in NS3/4A, NS5A, and NS5B. The role of this research was to screen for novel potential NS5B inhibitors from the cannabis compound database (CBD) using Deep Learning.

Methods:

Virtual screening of the CBD compounds was performed using a trained Graph Neural Network (GNN) deep learning model. Re-docking and conventional docking were used to validate the results for these ligands since some had rotatable bonds >10. About 31 of the top 67 hits from virtual screening and docking were selected after ADMET screening. To verify their candidacy, 6 random hits were taken for FEP/MD and Molecular Simulation Dynamics to confirm their candidacy.

Results:

The top 200 compounds from the deep learning virtual screening were selected, and the virtual screening results were validated by re-docking and conventional docking. The ADMET profiles were optimal for 31 hits. Simulated complexes indicate that these hits are likely inhibitors with suitable binding affinities and FEP energies. Phytil Diphosphate and glucaric acid were suggested as possible ligands against NS5B.

Untitled section

Keywords: Graph neural network, deep learning, Hepatitis C virus, Cannabis sativa, molecular simulation dynamics

Article notes

Untitled section

Received 2024 Sep 8; Accepted 2024 Nov 27; Collection date 2024.

Author notes

Introduction

The primary cause of chronic hepatitis, cirrhosis, and hepatocellular carcinoma (HCC) is the hepatitis C virus (HCV). Approximately 3% of the global population have a chronic HCV infection, and for 10 to 30 years, 30% of carriers are predicted to experience significant liver-related illnesses, such as HCC.¹ The number of cases of acute hepatitis C that have been reported has increased by 129% since 2014.²

HCV is an enveloped +ssRNA hepacivirus of the Flaviviridae family with 7 known genotypes and unknown genotypes with several subtypes that mostly affect the liver. Its 9.6 kb genome is flanked by 5′ and 3′ UTRs with a single ORF that encodes for a polyprotein with 3000 amino acids³ and is post-translationally processed by cellular and viral proteases to yield 11 viral proteins composed of 7 nonstructural proteins (NS) which are the 2 proteins necessary for the formation of the virion (p7 and NS2), as well as 5 proteins that make up the cytoplasmic viral replication complex (NS3, NS4A, NS4B, NS5A, and NS5B), and 3 structural proteins (1 nucleocapsid protein and 2 envelope proteins). Also encoded by the core region is an alternative open reading frame protein (ARFP) or F protein, whose function is still unknown⁴ even though some roles have been proposed including; regulation of Viral Translation, induction of apoptosis, interaction with host factors, and in viral morphogenesis.

NS5B is a 66 kDa heart-shaped catalytic component of the HCV replication complex, and owing to selectivity and non-toxicity, it makes a good therapeutic target due to the lack of mammalian counterparts. Like other polymerases, it has the palm, thumb, and finger domains surrounding the enzyme active site in the palm domain.⁵ Many anti-HCV Direct Acting Antivirals (DAAs) that are either nucleotide inhibitors (NIs) or non-nucleotide inhibitors (NNIs) targeting NS5B are in development.^6,7 At least 4 allosteric inhibitory sites have been reported (thumb site I or II, palm site I or II)⁸ with their mechanisms of action detailed elsewhere.^5,9
-13 However, these allosteric inhibitors have several drawbacks including low potency in enzymatic assays,¹⁴ lack of cellular potency during an HCV sub-genomic replicon assay,¹⁵ high lipophilic character, and low genetic barrier to resistance.¹⁵ It is possible for distinct forms of HCV with various amino acid changes that confer treatment resistance to coexist inside the same host in the context of quasispecies. While some alterations are not linked to drug resistance, others can result in a phenotypic decrease in susceptibility to 1 or more antiviral drugs.¹⁶ documented several substitutions with no discernible antiviral resistance to Dasabuvir and Sofosbuvir. Multi-drug-resistant RAASs variants of NS3/4A in GP1 and GP5 along with DAA-specific NS3/4A, NS5A, and NS5B were identified pan-genotypically.¹⁷ Resistance to daclatasvir and sofosbuvir was conferred by L2003M and S2702T of NS5B, respectively. D1194A NS3/4A was triple DAA (simeprevir, faldaprevir, and asunaprevir) resistant. Double-drug resistant variants included R1181K (faldaprevir and asunaprevir), A1182V and Q1106K/R (faldaprevir and simeprevir), T1080S (faldaprevir and telaprevir), while single drug-resistant variants were V1062L (telaprevir), D1194E/T (simeprevir), D1194G (asunaprevir), S1148A/G (simeprevir), and Q1106L (Boceprevir) of NS3/4A were determined.¹⁷ Many NS5B RAASs in genotypes 1, 3, 4, and 5 have been reported.¹⁷ Other mutations in NS5B that were related to DAAs resistance include E237G, S282R, L320F, V321A, and V321I.¹⁸ S282T induces high sofosbuvir-resistance, Q309R is a ribavirin-associated resistance, E237G was identified in the successfully amplified non-responder sample.¹⁹

With a multiverse of biochemical compounds (cannabinoids and non-cannabinoids (Phenolics, Terpenes, and Alkaloids)), cannabis spp. has been reported by several communities to have medicinal activities against several disease-causing pathogens²⁰ including HCV. HCV Population had lower rates of diabetes and obesity when they consumed cannabis, however, whether this translates into lower mortality should be investigated.²¹ It has been observed to address key challenges like nausea, and depression in HCV Patients. Some individuals receiving therapy for HCV may benefit virologically and symptomatically from modest cannabis use by helping them stick to the difficult medication regimen.^22,23 Faster decay of HIV RNA among cannabis users who were HCV co-infected and reduced risk of steatosis was observed.^24,25 The advantages of using cannabis for treating HCV from a biological and clinical standpoint, as well as the efficacy of this treatment, should be investigated in bigger study populations. Since there is no known vaccine or treatment for this virus, there is an urgent need for effective ways to manage and eventually eradicate this illness. This work aimed to investigate, using in-silico methodologies, if cannabis chemicals could block NS5B from HCV. Artificial intelligence (AI) was combined with traditional docking, binding energy, and simulation studies to validate AI results, quantify protein-ligand binding, and forecast the stability of the complexes. This combination significantly decreases the time and increases the accuracy needed for drug development.

Molecular dynamics (MD) simulations have become an invaluable tool in the fields of biology and drug development. These computational techniques allow researchers to study molecular systems at an atomic level over time, providing insights that are otherwise difficult to obtain through experimental approaches alone. By simulating the behavior of biomolecules, MD offers several key advantages in understanding biological processes and accelerating the drug discovery pipeline.^26
-28 Molecular dynamics (MD) simulations offer critical advantages in biology and drug development. They provide atomistic insights into the behavior of biomolecules, allowing researchers to observe dynamic processes like protein folding and conformational changes, which are difficult to capture experimentally.²⁹ In drug discovery, it helps predict ligand binding affinities, guiding the identification of potential drug candidates by simulating interactions between drugs and target proteins.³⁰ Additionally, MD is valuable for studying protein stability and misfolding, shedding light on diseases like Alzheimer’s.

MD simulations also allow for the investigation of complex biological processes, such as enzyme catalysis, and help researchers explore how proteins function under different conditions, like varying solvent environments. This flexibility enhances understanding of biological mechanisms and supports the design of drugs that better target disease processes. Furthermore, MD is cost-effective and accelerates research by reducing the need for exhaustive lab experiments. These benefits make MD simulations indispensable for advancing biological research and drug development.²⁷

Materials and Methods

Structures

The protein structure was downloaded from pdb (PDB ID: 3FQL) while cannabis small molecules were downloaded from the Cannabis Compound Database (CBD).

Geometric deep learning virtual screening

PDBbind database³¹ Protein-ligand complexes were used as inputs to train the model. Complexes that were also part of CASF-2016³² and those that failed pre-processing likely due to errors such as incomplete structure data, missing ligand coordinates, or errors in molecular surface generation, were excluded. The remaining complexes were randomly divided into training (14 000) and test (2367) sets. The detailed representation of ligand molecules is as in Méndez-Lucio et al.³³ and the protein targets were processed using a pipeline described by Gainza et al.³⁴ Two separate residual graph convolutional neural networks with the same architecture, 1 for the ligand and the other for the target, were used to extract features to build the model, the extracted features were concatenated and used to build a mixed-density network (MND).

Initially, a linear layer is used to project the node and edge features to a 128-dimensional embedding. Each node and edge was updated using a series of 3 GNNs depending on the nodes that were next to them and the kinds of edges that connected them. The GNN initially updates each edge in the graph by using a multi-layer perceptron (MLP) on the concatenation of the edge features and the features of the 2 connecting nodes and the updated edge features are used to update the node feature.³³ The modified edge and node features can be utilized as input for a subsequent convolution round because they contain information about not only the core atom but also its surrounding neighbors. Three convolutions were used. The node and edge features were then processed by the remaining GNN blocks.³³

After being pairwise concatenated, the node features recovered by the GNNs and residual GNNs for the ligand and target are fed into an MND.³⁵ Concatenated target and ligand node information are combined to construct a hidden representation by the MND using an MLP. The MND’s outputs are computed using the hidden representation. Furthermore, by connecting neighboring nodes, the retrieved ligand node properties were utilized to forecast the type of bond and atom, aiding in the learning of molecular structures and speeding up the training. Every MLP that is utilized consists of a linear layer, and an ELU activation function. The Exponential Linear Unit (ELU) was chosen over ReLU due to its benefits in stabilizing training by allowing for smoother gradients, especially in complex networks like GNNs. ELU’s ability to output negative values helps combat the vanishing gradient problem, which can be crucial in deep architectures. Experimenting with the RELU Section has been added to future work. A dropout rate of 0.1 was employed³³ appropriate to prevent over-regularization while maintaining model generalization. The dropout rate of 0.1 was determined based on previous studies and experimentation with similar deep-learning models. Given the complexity of GNNs and the molecular nature of the data, a smaller dropout rate was appropriate to prevent over-regularization while maintaining model generalization. Experimentation confirmed that 0.1 provided the right balance between reducing overfitting and maintaining predictive accuracy.

Training

The Adam optimizer was utilized to update the model weights at a learning rate of 0.002. The loss function was minimized during model training. 16 protein-ligand complexes were used as the batch size for 150 epochs of training the model. The learning rate (0.002), batch size (16), and number of epochs (150) were primarily based on best practices from the literature on geometric deep learning models, including those for graph neural networks (GNNs). These values have demonstrated stable convergence with similar contexts, and they balance training time and performance and were refined based on preliminary runs and validation accuracy. A potential specific to a given target-ligand complex was defined using the loss function. This potential was then used to score the target-ligand complex’s 3-dimensional structure by summing over all possible pairs, calculating the negative log-likelihood for each target-ligand node pair, and calculating the distances separating each target node from each ligand node in that particular conformation. This negative log-likelihood minimizes deviations between predicted and true binding conformations, making it ideal for the drug design domain where precision in ligand positioning is crucial. The likelihood of finding the target-ligand combination in that particular conformation increases with a decreasing value.^32,33

Benchmarking

The CASF-2016 benchmark,³² which includes 285 carefully chosen protein-ligand complexes, was used to evaluate this. The preprocessing of the structures from this benchmark was identical to that of the training set. The power of screening and scoring was assessed.³³

Prediction of binding conformations

The relative location of the ligand in Euclidean space, the dihedral angles of all rotatable bonds in the molecule, and the Euler angles were used to represent the ligand conformation as a vector. Differential evolution³⁶ was used to determine which ligand conformation would interact with the target binding site most likely following the model, that is, minimize the potential learned by the model for that particular complex. Using a population size of 150, the global optimization was performed up to 500 iterations with a recombination constant of 0.8 and a mutation constant randomly modified from 0.5 to 1.0 in every generation. Euler angles and the dihedrals of rotatable bonds were limited to values between −π and π without seeding.³³ The detailed geometric learning protocol can be accessed in the original publication.³³ The deep neural network learns the parameters of a mixture model that is employed as a probability density function. This probability density function is used to determine the most likely distance separating a ligand atom from a specific point in the molecular surface of the binding site. The potential is determined as the combination of the negative log-likelihood of all pairwise combinations of ligand atoms and points in the molecular surface. The optimal conformation is the one that minimizes the potential, that is, the ligand conformation in which every atom is separated from the target surface by the most likely distance.

Conventional Docking

Structure preparation

The protein structure was cleaned and preprocessed by assigning bond orders using the CCd database, adding Hydrogens, creating zero bond orders, and creating disulfide bonds as well as generation of het states using Epik³⁷ (PH 7 ± 2 units), water, and ligands were removed and H-bond assignment was done using PROPKA in Maestro 12.8 in the protein preparation wizard. The ligands were prepared by Ligprep.³⁸ Briefly, the ligands were imported into the maestro workspace, and only those with a maximum of 300 atoms were subjected to OPLS4 Force force field,³⁹ and ionization was done using Epik³⁷ to generate possible states at PH = 7±2 units.

Glide docking

The prepared structures were subjected to XP Glide⁴⁰ with only ligands of atoms and rotatable bonds equal to or below 300 and 100 respectively were selected with a Vander Waals scaling factor of 0.8 and a partial charge cutoff of 0.15. The receptor was rigid with flexible ligand sampling and to sample nitrogen inversions and ring conformations. Bias sampling torsion was set for all predefined functional groups and Epik state penalties were added to docking scores. Post-docking minimization was performed with 10 poses per ligand with a threshold of rejecting minimized pose set to 0.5 kcal/mol.

Binding affinity determination

SeeSAR 12.1.0⁴¹ was used to perform binding affinity calculations. A binding site was defined by the co-crystallized ligand in the receptor PDB file and copied to the docking mode after the generated protein and docking library were loaded into the Biosolveit workspace.⁴¹ Docking calculations were done for each compound in a standard docking mode with defaulted settings and parameters. The affinity of the generated poses was then assessed and the best poses were selected based on these affinities.

Binding energy is calculated from the HYDE score function (equation (1)), which relies on intrinsically balanced terms of atom-specific desolvation, hydration, and hydrogen bonds based on the logP atomic increment system,⁴² while also taking into account the Torsion angle values to the binding conformation of the protein-ligand. The quotient of G and the number of non-hydrogen atoms in the molecule is used to define a ligand’s binding affinity (equation (2))⁴³ and expressed as.

(1)

Δ GHyde = Σ atom i [\begin{matrix} Δ Gi Dehydn \\ + Δ Gi h - bonds \end{matrix}]

(2)

LE = Δ G ⁄ N

where ΔG = −RT ln Ki and N = number of non-hydrogen atoms.

The selection of the best poses was based on their visual HYDE scores while also considering a statistics-based torsional analysis.

ADMET screening

All compounds were subjected to SeeSAR,⁴¹ AdmetLab2.0,⁴⁴ and QikProp.^29,45 The protocol used for SeeSar is as in the Binding Energy determination above. The concatenated sdf file was uploaded to the AdmetLab2.0 website (https://admetmesh.scbdd.com/) in the ADMET screening mode. It employs a Multi-task Graph Attention (MGA) framework made up of input, Relation graph convolution network (RGCN) layers, attention layer, and fully connected (FC) layers.

Molecular simulation dynamics

All simulations were done in GROMACS 2021.4.⁴⁶ The protein topology was generated using AMBER99SB⁴⁷ force field from pdb2gmx module while ligands were parametrized by using Generalized Amber Force Field (GAFF2)⁴⁸ at antechamber website. The complexes were placed in the hexahedral box using gmx editconf, solvation was done using gmx solvate, TIP3P water molecules were added to the systems, and Na⁺ and Cl⁻ ions were added using genion to a concentration of 0.15 m. Equilibration was done with steps set to 20 000 in 10 ns at 298 K and 1 atm with position restraints Force Constant set to 700 kJ/mol. This allowed water molecules and ions to move freely. Production MD simulations lasted for 100 ns and they were monitored by checking system energies during the simulations. Pymol, vmd, and Gromacs binaries⁴⁶ were used for the analysis of the results.

MM/G(P)BSA calculations

Three energetic terms are taken into account in the computation of binding energy when determining the free energy of complex formation in conjunction with MD simulations: (1) variations in the system’s potential energy in a vacuum; (2) polar and apolar solvation of the various species; and (3) the entropy related to complex development during the gaseous phase.

The Uni-GBSA-based tool unigbsa-traj was utilized to do MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) calculations automatically for every simulated system.^45,49 The formula for estimation of free energy is explained in detail in Charles et al,⁴⁵ Wang et al,⁵⁰ and Gilson and Honig.⁵¹

Free energy perturbation (FEP) calculations

The CHARMM-GUI Free Energy Calculator and NAnoscale Molecular Dynamics (NAMD) were used to calculate absolute free energy.^29,52,53 The simulated ligands (docked postures) were uploaded as a single concatenated SDF file, and CHARMM General Force Field (CGenFF v1.x)⁵⁴ was used for ligand parametrization and topology construction to build the NAMD inputs and post-processing scripts. Counter ions (KCl) were used to neutralize all of the systems that were collected to produce input files and post-process scripts. Applying restraining potentials to limit the ligand’s location in a receptor during FEP/MD allowed for the calculation of binding free energy using the double decoupling approach. With the resulting inputs, the TIP3P water model and Langevin piston pressure were applied to the system’s NPT ensemble at 300 K and 1 atm. The FE values were captured in history files obtained via FEP lambda replica exchange MD (λ-REMD) using simple overlap sampling (SOS). By measuring the FE values throughout the last 6 ns, the final FE values from 10 ns FEP/λ-REMD simulations were computed.²⁹ The sequence of the used methods is shown in Figure 1.

Figure 1.The sequence of the steps followed in this study.

Results and Discussion

Geometric deep learning virtual screening

The results show that all optimal binding conformations had negative values (Figure 2 and Table 1) and only the best 100 molecules were considered for further studies based on their potential scores, with malic acid having the best score of −286.78, however, the complexes of those ligands with many rotatable bonds (⩽10) could not achieve successful optimizations and terminated after 500 iterations (Figure 2), this is because an increase in rotatable bonds is associated with the inefficiency of the optimization algorithms when working with a large number of degrees of freedom. These ligands justified the need for conventional docking methods integration, all compounds having rotatable bonds above 10 underwent the same re-docking process⁵⁵ and average potential scores were reported in Table 1.

Figure 2.The results of geometric deep learning showing the distribution of some key parameters. Fun: statistical potential, nit: number of iterations, rmsd: root mean square deviation.

Table 1.The summary of the results of virtual screening, molecular docking, and ADMET studies. Fun: A potential specific to a given target-ligand complex, SA Score: Synthetic Accessibility Score, LC50FM: 96-h fathead minnow LC50, FDAMDD: maximum recommended daily dose, Ames: Mutagenicity test, SR-p53: Probability of being p53 actives.

Generic name	Fun	Docking score (kcal/mol)	Glide ligand efficiency ln	Glide XP Hbond	QikProp stars	Binding affinity range (mM)	SA score	LC50FM	FDAMDD	Ames	SR-p53
Phytyl diphosphate^r	−197.87	−12.038	−2.811	−2.000	3	0.198-19.699	4.851	5.102	0.905	0.004	0.004
Apigenin-7-o-p-coumarylglucoside^***	−256.98	−12.038	−2.541	−4.427	5	1.2E6-1.2E9	4.081	5.818	0.122	0.686	0.969
Isocitric acid	−204.56	−11.050	−3.100	−2.480	3	1E3-1E5	3.356	3.261	0.006	0.021	0.004
LPA (18:0/0:0)^*^r	−187.45	−11.050	−2.530	−2.384	3	139-13 835	3.463	3.18	0.765	0.015	0.009
Glucaric acid	−263.768	−10.667	−2.931	−2.664	2	57-5628	4.092	2.556	0.004	0.02	0.004
LPA (18:1(11Z)/0:0)^*^r	−242.98	−10.621	−2.432	−2.102	3	70 986-7 052 884	3.664	2.805	0.887	0.024	0.008
LPA (18:2(9Z,12Z)/0:0)^*^r	−214.19	−10.210	−2.338	−1.310	3	0.768-76	3.846	3.58	0.947	0.004	0.01
Sativic acid^r	−194.50	−10.062	−2.408	−1.777	3	19-1933	3.524	3.041	0.019	0.028	0.017
LPA (18:1(9Z)/0:0)^r	−216.31	−9.978	−2.285	−2.568	3	1.515-150	3.664	3.257	0.889	0.006	0.009
Cynaroside^**^p	−196.45	−9.935	−2.225	−3.695	2	29-2868	3.924	4.845	0.031	0.757	0.835
Cannabisin A^p	−222.58	−9.873	−2.064	−2.558	5	34 041-3 382 134	2.817	4.116	0.493	0.494	0.809
Quercetin-o-glucoside^**^p	−277.72	−9.621	−2.140	−5.511	4	6179-613 939	4.008	4.785	0.02	0.809	0.781
Gluconic acid	−268.92	−9.602	−2.694	−3.600	2	3393-337 072	4.06	−0.392	0.002	0.053	0.003
Cannabitriol	−264.72	−8.908	−2.112	−1.184	0	2728-2 710 000	3.799	5.268	0.82	0.113	0.568
Malic acid	−286.78	−8.808	−2.755	−1.628	4	596-59 228	2.893	2.579	0.011	0.02	0.004
Quinic acid	−222.48	−8.736	−2.450	−3.300	2	15 396-1 529 705	3.598	0.807	0.016	0.031	0.003
cannabisin d^r	−222.08	−8.722	−1.806	−2.754	4	7398-735 004	3.412	5.418	0.946	0.791	0.811
Orientin^p	−256.86	−8.627	−1.932	−4.613	5	947-94 045	4.072	5.014	0.012	0.822	0.714
Isocannflavin B	−204.54	−8.540	−1.988	−1.567	0	3103-308 292	2.724	5.884	0.303	0.521	0.893
Apigenin-7-o-glucoside^**	−256.76	−8.522	−1.922	−1.790	1	83 708-8 316 933	3.804	4.804	0.019	0.693	0.871
N-trans-Caffeoyltyramine^p	−266.49	−8.499	−2.077	−1.912	0	8.240-819	2.061	4.611	0.64	0.386	0.82
3-[2-(3-Isoprenyl-4-hydroxy-5-methoxy-phenyl)ethyl]-5-methoxyphenol	−209.68	−8.456	−2.004	−1.296	1	269-26 722	2.454	5.3	0.827	0.017	0.845
Cannabistilbene I	−196.04	−8.378	−2.026	−0.136	0	5.135-510	2.361	5.545	0.917	0.023	0.465
Quercetin^p	−201.67	−8.360	−2.043	−2.400	0	9.593-953	2.545	5.222	0.31	0.657	0.888
Cannabidiolic acid	−228.47	−8.282	−1.945	−1.180	0	185-18 368	3.626	5.699	0.787	0.023	0.513
Arachidic acid^r	−148.87	−8.237	−2.013	−1.158	5	10.4-1037	1.665	4.743	0.014	0.005	0.034
Etofenoprox	−176.93	−8.220	−1.897	0.000	2	0.440-44	2.025	6.485	0.782	0.186	0.005
Delta-9-tetrahydro cannabinolic acid A	−279.84	−8.202	−1.926	−1.380	1	1028-102 145	3.617	5.55	0.689	0.012	0.413
Resmethrin^*	−232.34	−8.053	−1.909	0.000	1	11-1133	3.317	7.826	0.631	0.022	0.002
Astragalin^**	−222.78	−8.016	−1.795	−2.564	2	3225-320 374	3.884	4.732	0.011	0.775	0.853
Isovitexin	−197.23	−8.002	−1.805	−2.446	3	619 407-61 541 783	4.081	5.818	0.122	0.686	0.969
Sofosbuvir^r	−201.60	−6.887	−1.503	−0.494	0	64-6368	4.375	4.036	0.9	0.371	0.02

Conventional docking and binding affinity determination

The results of docking and binding affinity determination are shown in Table 1 and they agree with the results of geometric learning, validating the compounds under investigation. Phytyl diphosphate had the best score of −12.275 kcal/mol, but iso-citric acid had the best ligand efficiency of −3.100. The compounds had acceptable values of Binding affinities⁴¹ with Isovitexin having the highest value of 619 407 to 61 541 783 nM (Table 1).

ADMET screening

The results for some important selected properties are presented in Table 1, and the summary of all properties is shown in Figure 3 and Supplemental File 2, because the presence of PAINS alert is not enough to justify the elimination of a hit candidate,⁴⁵ Compounds with pain alerts were only eliminated if at least other 4 ADMET properties were out of range (except the number of bonds, which was not considered as an ADMET property). The Synthetic Accessibility Score (SA Score) quantifies how simple it is to synthesize drug-like compounds. A score of less than 6 indicates that the compound is simple to synthesize.^29,44 An estimate of the hazardous dosage threshold of substances in humans can be obtained from the maximum recommended daily dose (FDAMDD). The likelihood of being toxic is the output value, and it ranges from 0 to 0.3 for excellent (less likely to be harmful); 0.3 to 0.7 for medium; and 0.7 to 1.0 for poor (more likely to be toxic).⁴⁴ The mutagenicity test is the Ames test. The most used test for determining a compound’s mutagenicity is the mutagenic effect, which closely correlates with carcinogenicity. The numbers indicate the probability of being harmful; 0 to .3 indicates excellent (less likely to be mutagenic), .3 to .7 indicates medium, and .7 to 1.0 indicates poor (more likely mutagenic).⁴⁴ A fathead minnow LC50FM is defined as the concentration of the test chemical in water which results in 50% of the minnows dying after 96 hours. For LC50FM, the unit is −log10[(mg/L)/(1000 × MW)].⁴⁴ One reliable sign of DNA damage and other cellular stressors is the activation of p53. The output values of SR-p53 are the probabilities of being active, with 0 to 0.3 representing excellent (likely inactive), 0.3 to 0.7 representing medium, and 0.7 to 1.0 representing poor (most likely active).⁴⁴ The rest of the computed ADMET properties are presented in Supplementary File 2^41,44 and their tSNE distribution is in Figure S1.

Figure 3.The K-means cluster analysis of ADMET properties.

Molecular dynamics simulation studies

The potential of 6 protein-ligand complexes as NS5B inhibitors was demonstrated by the formation of stable complexes by each of the simulated complexes. The analysis of Root Mean Square Deviation (RMSD) evaluates the long-term structural stability of biomolecular simulations. It estimates the typical difference between an atom’s location in a molecular structure and its reference structure.^29,56 The RMSDs of each system show that during the experiment, every simulated complex reached stability (Figure 4A and S2A). Comparatively, the glucaric acid complex was more stable than the other compounds. The observed RMSD differences between the complexes and their corresponding proteins exhibit a variation that is driven by the ligand. When utilizing 1D RMSD, it is easy to believe that 2 structures that have the same RMSD from a reference frame are similar, but in practice, they can differ significantly. Alternatively, significantly more information can be obtained by computing the RMSD of each frame in the trajectory to all other frames in the other trajectory to give 2D RMSD.⁴⁵ Pairwise RMSDs of each trajectory were calculated to itself and the results are shown in Figure S5, the diagonal represents 0 (RMSD of a structure to its self), and all the structures had RMSDs ⩽ 1.8 Å, over simulation showing that the complexes were stable over simulation and had nearly same states but not identical except Phytyl diphosphate complex had more re-visited states than other complexes (Figure S5). These results agree with 1D RMSD results and show that the glucaric acid complex was more stable evidenced by lower RMSD values.

Figure 4.The molecular dynamics simulation results: (A) illustrates the RMSD variation for proteins across different complexes, (B) displays the total interaction energy for each system throughout the simulation, (C) shows the RMSF variation over time, and (D) presents the changes in the radius of gyration (Rg) during the simulation.

Root Mean Square Fluctuation (RMSF) analysis is another way to understand the flexibility and dynamic behavior of individual atoms or residues within a biomolecular system as well as their contribution to the flexibility of the whole molecule.^26,30,57 The main interacting residues have minimum RMSF values during simulation time compared to the ligand-free protein, supporting their stability and interactions with the simulated compounds, while the RMSF of non-interacting residues shows somewhat greater oscillations (Figure 4C). The distribution of the rmsf further supports the stability of all complexes (Figure S4D).

The Radius of Gyration (Rg) is a measure of the compactness or spread of a biomolecular structure in 3-dimensional space is valuable for analyzing the overall shape and structural fluctuations of biomolecules.^30,45 A relatively constant Rg value fundamentally signifies a stably folded structure and a reduction in Rg signifies an increase in stability.^29,30,45 All of the simulated systems had Rg values between 23.7 and 24.2 Å, indicating their stability. The simulation revealed a fairly progressive reduction in the gyration radius over time, indicating a gradual increase in the compactness and stability of all systems (Figure 4D). The distribution of Rg also shows that the Glucaric acid complex is more stable relative to other complexes (Figure S3C).

Solvent Accessible Surface Area (SASA) is a measure of the surface area of a biomolecule that is accessible to solvent molecules.⁵⁸ It plays a crucial role in analyzing the interactions between biomolecules and their surrounding solvent environment.⁴⁵ A higher value suggests an increase in the protein’s volume, indicating lower stability, whereas stable proteins typically exhibit minimal fluctuation throughout the simulation. The binding of a small molecule can alter the solvent-accessible surface area (SASA) and significantly impact the protein’s structure.^28,45,59 The SASA consistently decreased across all systems during the 100 ns simulation (Figure S2B). This decrease in SASA indicates a rise in compactness, consequently indicating enhanced stability across all systems. The parallel patterns observed in both SASA and Rg affirm the validity of the molecular dynamics simulation outcomes.⁴⁵ Their distribution also shows that all systems were stable over simulation with the glucaric acid complex being more stable than others (Figure S3C and S3D).

Non-bonded Molecular Mechanics (MM) interaction energy between the Ligands and their receptor was calculated to evaluate the magnitude of the interaction between the ligand and the protein.^60,61 The total interaction energy for all the systems was negative over 100 ns simulation with the complex of phytyl-diphosphate having the lowest energy values (Figure 4B and S3B). Vander Waals forces are short-range interactions that include London dispersion forces and dipole-dipole interactions, they enable the close interaction between the nonpolar regions of the ligand and the protein, ensure complex shape complementarity, predict the strength of the protein-ligand interaction with higher energy indicating stronger binding.⁶² Apigenin coumaryl glucoside had the best values while that of Isocitric acid was the worst, a fact justified by their size and proximity to the protein groups (Figure S9). The stability of the complex and total binding energy can be enhanced by the long-range electrostatic interactions, which can direct protein molecules toward their pre-binding orientations.⁶³ Phytyl Diphosphate had the best values of electrostatic energy while apigenin coumaryl glucoside had the worst values since the former has phosphate groups which results in strong electrostatic interaction energy. These energies are good predictions of binding affinity since they are considered while computing those affinities. The simulated complexes remained stable throughout the simulation according to this data (Figure S2C and S2D). But it’s crucial to remember that this quantity is neither a binding energy nor a free energy.⁶⁰ Figure 5 displays the specifics of the molecular interactions that occur in the middle of the simulation. The summary and distribution of all these energy terms are shown in Figures S3B, S4A, S4B, and S4C.

Figure 5.An overview of the interactions between the selected compounds and the target protein. The interacting residues and compounds are shown in licorice representation, with the residues highlighted in green and the compounds in pink, while the rest of the protein is displayed in a cartoon format.

During the simulation, the minimum distances between active site residues and ligands were calculated. Since most of the time, these minimum distances were less than 3 Å, a conventional hydrogen bond can still form as long as the acceptor and donor are orientated correctly (Figure S3A).

PCA (Principal Component Analysis)^27,64,65 was performed particularly on alpha carbon data from the last 25 ns to acquire insights into the dynamic behavior of both the complexes, incorporating structural and energy data. This analysis aimed to explore the conformational range of the complexes, distinguishing various regions within the energy landscape explored during the MD simulation. The complexes defined discrete conformational clusters and suggested stability by occupying compact subspaces. A graph illustrating the motion and displacement of atomic fluctuations within the complexes was created by utilizing eigenvectors 1 and 2. The first few eigenvalues had greater values, and the remaining eigenvalues were in a declining order.⁴⁵ All complexes showed PC ranging from ~15 to ~−15 (Figure S6). The simulation results of the apo-protein are shown in Figure S10. Interacting residues have increased positive correlation, as indicated by the Dynamic Cross-Correlation (DCC) data, indicating that they are interacting with the target. Residues located in the active site showed slightly elevated positive cross-correlation peaks upon ligand binding, suggesting a high occupant binding affinity (Figure S8).⁴⁵

The average binding free energy of simulated hits was estimated using MM/P/GBSA calculations, and the outcomes are in good agreement with other studies. Since all energies were negative, proteins and ligands were strongly bound together (Table 2). Additionally, FEP/MD simulations were performed; Table 2’s results show that these ligands have sufficient binding affinities, making them eligible for in-vitro research. To evaluate the convergence and dependability of the findings, the FE values were averaged using the standard error of the mean (SOS). The SEM of all data is less than 0.5 kcal/mol, indicating that all systems have converged in less than 10 ns.^29,52,66

Table 2.The table shows the results of MM_G/PBSA and FEP/MD energy calculations.

Ligand	MM_GBSA (kcal/mol)	MM_PBSA (kcal/mol)	ΔG FE (kcal/mol)
Phytyl Diphosphate	−57.331	−31.784	−9.45
Apigenin Coumarylglucoside	−46.342	−20.304	−8.24
Isocitric Acid	−43.111	−29.160	−9.06
LPA (18:0/0:0)	−60.753	−32.981	−8.64
Glucaric Acid	−66.876	−35.834	−7.32
LPA (18:1(11Z)/0:0)	−53.496	−28.507	−7.86

Because of many RAASs in HCV proteins, reported DAAs resistance-enhancing RAASs were introduced in this protein with several combinations as shown in Table S1 and these mutants were re-docked with the identified compounds. 7 ligands still had significant docking scores with Glucaric acid having the highest values (Figure 6 and Table S2). Paired two-tailed t-tests to compare the wild type and mutants are also shown in Tables S4, S5, S6, and S7. These results show that the docking scores of many ligands are significantly reduced except for the 7 ligands. This may suggest that some of these compounds, with glucaric acid leading, may still bind the resistant phenotypes. While these docking results provide valuable insights, they alone are insufficient to fully demonstrate the efficacy of the identified ligands against the mutant proteins. Additional investigations, such as molecular dynamics simulations of the mutant complexes, will be necessary to further validate and understand their interactions and potential effectiveness.^26,57,64,67

Figure 6.The results of docking the ligands into the active sites of mutants.

Studies have demonstrated the hepatoprotective, anti-inflammatory, cholesterol-lowering, anti-oxidant, and anti-carcinogenic properties of glucaric acid and its derivative D-saccharic acid-1,4-lactone.⁶⁸

With the pan genomic approach to HCV therapy, the candidates were also docked with NS3/4 protease (PDB ID: 3P8N), and the results are presented in Table S3, 4 ligands including Glucaric acid had docking scores better than −7.00 kcal/mol. Having inhibitors inhibiting both proteins will be a bonus for drug development. To further discover the global therapeutic world of glucaric acid, insight into the compound’s potential pharmacological targets and associated biological pathways was determined. This was done using Swiss Target Prediction to predict other possible targets of Glucaric Acid. Results (Figure 7) show that this compound may be active against several targets including the Neuronal acetylcholine receptor protein alpha-7 subunit (Probability = .765) and Squalene synthetase (Probability = .103). Inhibiting the former means that this compound can also be a competitive antagonist at a neuromuscular junction and hence can competitively compete with bungarotoxin minimizing its blockade activities at the junction⁶⁹ while inhibiting the latter is a potential therapeutic strategy for lowering cholesterol levels in individuals with hypercholesterolemia, as well as Anti-Cancer Potential Since cholesterol, is essential for the formation of lipid rafts and cell signaling pathways involved in cancer cell proliferation and survival, inhibiting squalene synthase can potentially inhibit tumor growth and metastasis.^68,70 The identified hits are shown in Figure S11.

Figure 7.The results show other possible targets for glucaric acid.

Some studies have also highlighted cannabis-delivered phytochemicals to inhibit this protein, 1 study highlights the use of terpenes, naturally occurring compounds found in plants including cannabis, as potential HCV NS5B polymerase inhibitors. The study identified several terpenes, such as mulberrofuran G, cochlearine A, and pawhuskin B, that demonstrated stable binding interactions with the NS5B enzyme. These terpenes were suggested as promising candidates for inhibiting viral replication, showing better docking results compared to traditional HCV drugs like sofosbuvir ,⁷¹ based on the distribution of molecular simulation parameters of complexes, our complexes show more stability. Additionally, another study explored various phytochemicals, including gallic acid, catechin, resveratrol, apigenin, and silibinin, which have antiviral properties. These compounds were docked against the HCV NS5B enzyme, with silibinin showing the best binding affinity, suggesting it could be a potent antiviral agent against HCV ,⁷² the compounds evaluated in this study exhibit superior docking scores and binding affinity ranges compared to compounds above from.⁷² These results emphasize the strong potential of the studied compounds for effective target interaction, with simulations further supporting their viability, unlike the lower-performing compounds from previous studies, which lack binding affinity data and simulated validation (Table 3).

Table 3.Comparison of docking scores, binding affinity ranges, and simulation data between compounds from this study and those from Shakya.⁷² The compounds from this study demonstrate stronger docking scores, broader binding affinity ranges, and simulation validation, highlighting their superior potential for target interaction.

Compound	Docking score (kcal/mol)	Binding affinity range (mM)	Simulation conducted	Source
Phytyl diphosphate	−12.275	0.198-19.699	Yes	This study
Apigenin-7-o-p-coumarylglucoside	−12.275	1.2 × 10⁶ to 1.2 × 10⁹	Yes	This study
Isocitric acid	−11.050	1 × 10³ to 1 × 10⁵	Yes	This study
Silibinin	−10.57	No data reported	No	Shakya⁷²
Apigenin	−8.75	No data reported	No	Shakya⁷²
Resveratrol	−8.14	No data reported	No	Shakya⁷²

Future work and limitations

Molecular dynamics simulations will be essential to further investigate the impact of NS5B mutations at the molecular level, as well as studying these proteins under different solvent conditions, such as urea. While the docking results provide useful preliminary insights, they are insufficient to confirm the efficacy of the identified ligands against mutant proteins.^65,73,74 Additional studies, including molecular dynamics simulations of the mutant complexes, will be necessary to validate these findings and better understand ligand-protein interactions for potential therapeutic applications.²⁶ Molecular dynamics (MD) simulations and in silico methods, including docking and AI models, offer powerful tools for studying biological systems and drug discovery, but they come with limitations. MD simulations are computationally expensive, especially for large systems or long timescales, requiring significant resources. The accuracy of these simulations depends on force fields, which may not always represent molecular interactions accurately. Additionally, simulations often oversimplify complex biological environments, leading to less realistic results. Moreover, in silico findings require experimental validation to ensure their relevance. These methods are most effective when used alongside laboratory experiments to complement and verify computational predictions.^75,76 More hyperparameter tuning as well as experimenting with ReLU will be a future optimization step to see if it improves the model’s convergence.

Conclusion

We have used geometric deep learning to screen for possible HCV inhibitors of NS5B from cannabis sativa natural compounds. These compounds’ data show that they may be potential DAAs against HCV. It is however very important to keep the RAASs in check, given that they represent the biggest challenge to HCV treatments, and screening them and recommending treatment post their identification will be a step ahead. The best hits identified were Glucaric acid and phytyl diphosphate, with the former retaining the ability to interact with mutants studied as well NS3/4 hence a possibility of polypharmacology and better than the used standard of sofosbuvir.

This work broadens the scope of cannabis research by shedding light on the potential therapeutic applications of the plant beyond its well-known effects, improving our understanding of its pharmacological properties, and highlighting its potential as a source of novel antiviral compounds.

This study further guides in vitro and in vivo experimental studies for validation of the actual efficacy, safety, and inhibitory mechanisms of these compounds. Cannabis-derived compounds may offer advantages such as natural sourcing, potentially fewer side effects, and diverse chemical structures for drug optimization and could also inform the design of more potent and selective HCV inhibitors.

Supplemental Material

Footnotes

Footnote Group

References

Untitled section

References

1. Shepard CW, Finelli L, Alter MJ. Global epidemiology of hepatitis C virus infection. Lancet Infect Dis. 2005;5:558-567.
2. CDC. 2021. Hepatitis C and Viral Hepatitis Surveillance Report. n.d. Accessed January 14, 2024. https://www.cdc.gov/hepatitis/statistics/2021surveillance/hepatitis-c.htm
3. Venkatesan A, Febin Prabhu Dass J. Review on chemogenomic approaches towards hepatitis C viral targets. J Cell Biochem. 2019;120:12167-12181.
4. Xu Z, Choi J, Yen TS, et al. Synthesis of a novel hepatitis C virus protein by ribosomal frameshift. EMBO J. 2001;20:3840-3848.
5. Hang JQ, Yang Y, Harris SF, et al. Slow binding inhibition and mechanism of resistance of non-nucleoside polymerase inhibitors of hepatitis C virus. J Biol Chem. 2009;284:15517-15529.
6. Moradpour D, Penin F. Hepatitis C virus proteins: from structure to function. Curr Top Microbiol Immunol. 2013;369:113-142.
7. Irekeola AA, Ear ENS, Mohd Amin NAZ, Mustaffa N, Shueb RH. Antivirals against HCV infection: the story thus far. J Infect Dev Ctries. 2022;16:231-243.
8. Sorbo MC, Cento V, Di Maio VC, et al. Hepatitis C virus drug resistance associated substitutions and their clinical relevance: update 2018. Drug Resist Updat. 2018;37:17-39.
9. Winquist J, Abdurakhmanov E, Baraznenok V, et al. Resolution of the interaction mechanisms and characteristics of non-nucleoside inhibitors of hepatitis C virus polymerase. Antiviral Res. 2013;97:356-368.
10. Yi G, Deval J, Fan B, et al. Biochemical Study of the comparative inhibition of hepatitis C virus RNA polymerase by VX-222 and filibuvir. Antimicrob Agents Chemother. 2012;56:830-837.
11. Caillet-Saguy C, Simister PC, Bressanelli S. An objective assessment of conformational variability in complexes of hepatitis C virus polymerase with non-nucleoside inhibitors. J Mol Biol. 2011;414:370-384.
12. Boyce SE, Tirunagari N, Niedziela-Majka A, et al. Structural and regulatory elements of HCV NS5B polymerase–β-loop and C-terminal tail–are required for activity of allosteric thumb site II inhibitors. PLoS One. 2014;9:e84808.
13. Han D, Wang H, Wujieti B, et al. Insight into the drug resistance mechanisms of GS-9669 caused by mutations of HCV NS5B polymerase via molecular simulation. Comput Struct Biotechnol J. 2021;19:2761-2774.
14. Stammers TA, Coulombe R, Rancourt J, et al. Discovery of a novel series of non-nucleoside thumb pocket 2 HCV NS5B polymerase inhibitors. Bioorg Med Chem Lett. 2013;23:2585-2589.
15. Eltahla AA, Luciani F, White PA, Lloyd AR, Bull RA. Inhibitors of the hepatitis C virus polymerase; mode of action and resistance. Viruses. 2015;7:5206-5224.
16. Hasanshahi Z, Hashempour A, Ghasabi F, et al. First report on molecular docking analysis and drug resistance substitutions to approved HCV NS5A and NS5B inhibitors amongst Iranian patients. BMC Gastroenterol. 2021;21:443.
17. Izhari MA. Molecular mechanisms of resistance to direct-acting antiviral (DAA) drugs for the treatment of hepatitis C virus infections. Diagnostics. 2023;13:3102.
18. Nguyen TK, Van Le D. Identification of NS5B resistance-associated mutations in hepatitis C virus circulating in treatment naïve Vietnamese patients. Infect Drug Resist. 2022;15:1547-1554.
19. Ahmed HR, Waly NGFM, Abd El-Baky RM, et al. Distribution of naturally-occurring NS5B resistance-associated substitutions in Egyptian patients with chronic hepatitis C. PLoS One. 2021;16:e0249770.
20. Ssemuyiga C. In-vitro determination of antibacterial activity of Cannabis Sativa against staphylococcus aureusIn-vitro determination of antibacterial activity of Cannabis Sativa against staphylococcus aureus. 2019. Accessed November 12, 2023. http://dissertations.mak.ac.ug/handle/20.500.12281/7027
21. Barré T, Bourlière M, Ramier C, et al. Cannabis use is inversely associated with metabolic disorders in hepatitis C-infected patients (ANRS CO22 hepather cohort). J Clin Med. 2022;11:6135.
22. Fischer B, Reimer J, Firestone M, et al. Treatment for hepatitis C virus and cannabis use in illicit drug user patients: implications and questions. Eur J Gastroenterol Hepatol. 2006;18:1039-1042.
23. Sylvestre DL, Clements BJ, Malibu Y. Cannabis use improves retention and virological outcomes in patients treated for hepatitis C. Eur J Gastroenterol Hepatol. 2006;18:1057-1063.
24. Marcellin F, Miailhes P, Santos M, et al. Cannabis use and plasma human immunodeficiency virus (HIV) RNA levels in patients coinfected with HIV and hepatitis C virus receiving antiretroviral therapy: data from the ANRS CO13 HEPAVIH cohort. Clin Infect Dis. 2020;71:2536-2538.
25. Nordmann S, Vilotitch A, Roux P, et al. Daily cannabis and reduced risk of steatosis in human immunodeficiency virus and hepatitis C virus-co-infected patients (ANRS CO13-HEPAVIH). J Viral Hepat. 2018;25:171-179.
26. Khan FI, Hassan F, Anwer R, Juan F, Lai D. Comparative analysis of bacteriophytochrome agp2 and its engineered photoactivatable nir fluorescent proteins pairfp1 and pairfp2. Biomolecules. 2020;10:1-22.
27. Khan FI, Wei DQ, Gu KR, Hassan MI, Tabrez S. Current updates on computer aided protein modeling and designing. Int J Biol Macromol. 2016;85:48-62.
28. Khan FI, Ali S, Chen W, et al. High-resolution MD simulation studies to get mechanistic insights into the urea-induced denaturation of human sphingosine kinase 1. Curr Top Med Chem. 2021;21:2839-2850.
29. Charles S, Mahapatra RK. Artificial intelligence based de-novo design for novel plasmodium falciparum plasmepsin (PM) X inhibitors. J Biomol Struct Dyn. 2023;1-16. doi: 10.1080/07391102.2023.2279700
30. Charles S, Edgar MP, Mahapatra RK. Artificial intelligence based virtual screening study for competitive and allosteric inhibitors of the SARS-CoV-2 main protease. J Biomol Struct Dyn. 2023;41:15286-15304.
31. Liu Z, Li Y, Han L, et al. PDB-wide collection of binding data: current status of the PDBbind database. Bioinformatics. 2015;31:405-412.
32. Su M, Yang Q, Du Y, et al. Comparative assessment of scoring functions: the CASF-2016 update. J Chem Inf Model. 2019;59:895-913.
33. Méndez-Lucio O, Ahmad M, del Rio-Chanona EA, Wegner JK. A geometric deep learning approach to predict binding conformations of bioactive molecules. Nat Mach Intell. 2021;3:1033-1039.
34. Gainza P, Sverrisson F, Monti F, et al. Deciphering interaction fingerprints from protein molecular surfaces using geometric deep learning. Nat Methods. 2020;17:184-192.
35. Bishop CM. Mixture density networks. 1994. Accessed January 4, 2024. https://www.semanticscholar.org/paper/Mixture-density-networks-Bishop/4cf3569e045993dfe090749f26a55a768684ab86
36. Storn R, Price K. Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces. J Glob Optim. 1997;11:341-359.
37. Shelley JC, Cholleti A, Frye LL, et al. Epik: a software program for pk(a) prediction and protonation state generation for drug-like molecules. J Comput Aided Mol Des. 2007;21:681-691.
38. LigPrep User Manual. Schrödinger press LigPrep 3.4 user manual. 2015.
39. Lu C, Wu C, Ghoreishi D, et al. OPLS4: Improving force field accuracy on challenging regimes of chemical space. J Chem Theory Comput. 2021;17:4291-4300.
40. Friesner RA, Murphy RB, Repasky MP, et al. Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein−ligand complexes. J Med Chem. 2006;49:6177-6196.
41. SeeSAR version 13.0.0; BioSolveIT GmbH, Sankt Augustin, Germany, 2023, www.biosolveit.de/SeeSAR
42. Brethon A, Chantalat L, Christin O, et al. New caspase-1 inhibitor by scaffold hopping into bio-inspired 3D-fragment space. Bioorg Med Chem Lett. 2017;27:5373-5377.
43. Schärfer C, Schulz-Gasch T, Ehrlich HC, et al. Torsion angle preferences in druglike chemical space: a comprehensive guide. J Med Chem. 2013;56:2016-2028.
44. Xiong G, Wu Z, Yi J, et al. ADMETlab 2.0: an integrated online platform for accurate and comprehensive predictions of ADMET properties. Nucleic Acids Res. 2021;49:W5-W14.
45. Charles S, Edgar MP, Kasoma NA. The hunt for antipox compounds against monkeypox virus thymidylate kinase and scaffolding protein leveraging pharmacophore modeling, molecular docking, ADMET studies and molecular dynamics simulation studies. Virol Mycol. 2023;12:1-14.
46. Bauer P, Hess B, Lindahl E. GROMACS 2022.1 manual. 2022. Accessed February 2, 2024. 10.5281/ZENODO.6451567
47. Lindorff-Larsen K, Piana S, Palmo K, et al. Improved side-chain torsion potentials for the amber ff99SB protein force field. Proteins. 2010;78:1950-1958.
48. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. Development and testing of a general amber force field. J Comput Chem. 2004;25:1157-1174.
49. Yang M, Wang D, Zheng H. UNI-GBSA: an automatic workflow to perform MM/GB(PB)SA calculations for virtual screening. n.d. Accessed Fbruary 4, 2024. https://hermite.dp.tech/
50. Wang J, Hou T, Xu X. Recent advances in free energy calculations with a combination of molecular mechanics and continuum models. Curr Comput Aided-Drug Des. 2006;2:287-306.
51. Gilson MK, Honig B. Calculation of the total electrostatic energy of a macromolecular system: solvation energies, binding energies, and conformational analysis. Proteins. 1988;4:7-18.
52. Zhang H, Kim S, Giese TJ, et al. CHARMM-GUI free energy calculator for practical ligand binding free energy simulations with AMBER. J Chem Inf Model. 2021;61:4145-4151.
53. Gumbart J, Hénin J, Chipot C. In silico alchemy: a tutorial for alchemical free-energy perturbation calculations with NAMD. 2023. https://www.semanticscholar.org/paper/In-silico-alchemy%3A-A-tutorial-for-alchemical-with-Harrison/8fbb201c76e1e06686bc54634d5b3a100955db8d
54. Soteras Gutiérrez I, Lin FY, Vanommeslaeghe K, et al. Parametrization of halogen bonds in the CHARMM general force field: improved treatment of ligand-protein interactions. Bioorg Med Chem. 2016;24:4812-4825.
55. Li H, Leung KS, Ballester PJ, Wong MH. Istar: a web platform for large-scale protein-ligand docking. PLoS One. 2014;9:e85678.
56. Martínez L. Automatic identification of mobile and rigid substructures in molecular dynamics simulations and fractional structural fluctuation analysis. PLoS One. 2015;10:e0119264.
57. Qausain S, Khan FI, Khan MKA. Conserved acidic second shell residue modulates the structure, stability and activity of non-seleno human peroxiredoxin 6. Int J Biol Macromol. 2023;242:124796.
58. Wang S, Khan FI. Investigation of molecular interactions mechanism of pembrolizumab and PD-1. Int J Mol Sci. 2023;24:10684.
59. Stephens DE, Khan FI, Singh P, et al. Creation of thermostable and alkaline stable xylanase variants by DNA shuffling. J Biotechnol. 2014;187:139-146.
60. Arantes PR, Polêto MD, Pedebos C, Ligabue-Braun R. Making it rain: cloud-based molecular simulations for everyone. J Chem Inf Model. 2021;61(10):4852-4856.
61. Charles S, Pius Edgar M, Mahapatra RK. Host-Directed Anti-Fusion Aptamers and Small Molecules as Respiratory Syncytial Virus (RSV) Inhibitors: An in silico-based study. J Biomed Biotechnol. 2024;7:485–497. Accessed February 10, 2024. 10.26502/jbb.2642-91280171
62. Bitencourt-Ferreira G, Veit-Acosta M, de Azevedo WF. Van der waals potential in protein complexes. Methods Mol Biol. 2019;2053:79-91.
63. Zhang Z, Witham S, Alexov E. On the role of electrostatics in protein-protein interactions. Phys Biol. 2011;8:035001.
64. Khan FI, Shahbaaz M, Bisetty K, et al. Large scale analysis of the mutational landscape in β-glucuronidase: a major player of mucopolysaccharidosis type VII. Gene. 2016;576:36-44.
65. Khan FI, Gupta P, Roy S, et al. Mechanistic insights into the urea-induced denaturation of human sphingosine kinase 1. Int J Biol Macromol. 2020;161:1496-1505.
66. Jiang W, Roux B. Free energy perturbation Hamiltonian replica-exchange molecular dynamics (FEP/H-REMD) for absolute ligand binding free energy calculations. J Chem Theory Comput. 2010;6:2559-2565.
67. Hassan F, Khan FI, Song H, Lai D, Juan F. Effects of reverse genetic mutations on the spectral and photochemical behavior of a photoactivatable fluorescent protein PAiRFP1. Spectrochim Acta. 2020;228:117807.
68. Ayyadurai VAS, Deonikar P, Fields C. Mechanistic understanding of D-glucaric acid to support liver detoxification essential to muscle health using a computational systems biology approach. Nutrients. 2023;15. doi: 10.3390/nu15030733
69. Chang CC, Chen TF, Chuang ST. Influence of chronic neostigmine treatment on the number of acetylcholine receptors and the release of acetylcholine from the rat diaphragm. J Physiol. 1973;230:613-618.
70. Seiki S, Frishman WH. Pharmacologic inhibition of squalene synthase and other downstream enzymes of the cholesterol synthesis pathway: a new therapeutic approach to treatment of hypercholesterolemia. Cardiol Rev. 2009;17:70-76.
71. Karpiński TM, Ożarowski M, Silva PJ, et al. Discovery of terpenes as novel HCV NS5B polymerase inhibitors via molecular docking. Pathogens. 2023;12:842.
72. Shakya AK. Natural phytochemicals: potential anti-HCV targets in silico approach. J Appl Pharm Sci. 2019;9:94-100.
73. Qausain S, Khan FI, Lai D, et al. Mechanistic insights into the urea-induced denaturation of a non-seleno thiol specific antioxidant human peroxiredoxin 6. Int J Biol Macromol. 2020;161:1171-1180.
74. Khan FI, Aamir M, Wei DQ, Ahmad F, Hassan MI. Molecular mechanism of ras-related protein Rab-5a and effect of mutations in the catalytically active phosphate-binding loop. J Biomol Struct Dyn. 2017;35:105-118.
75. Palsson B. The challenges of in silico biology. Nat Biotechnol. 2000;18(11):1147–1150. Accessed February 20, 2024. 10.1038/81125
76. Filipe HAL, Loura LMS. Molecular dynamics simulations: advances and applications. Molecules. 2022;27:2105.

Studies

Conditions

Cannabinoids and active compounds

Related research

Abstract

Introduction:

Methods:

Results:

Untitled section

Article notes

Untitled section

Author notes

Introduction

Materials and Methods

Structures

Geometric deep learning virtual screening

Training

Benchmarking

Prediction of binding conformations

Conventional Docking

Structure preparation

Glide docking

Binding affinity determination

ADMET screening

Molecular simulation dynamics

MM/G(P)BSA calculations

Free energy perturbation (FEP) calculations

Results and Discussion

Geometric deep learning virtual screening

Conventional docking and binding affinity determination

ADMET screening

Molecular dynamics simulation studies

Future work and limitations

Conclusion

Supplemental Material

Footnotes

Footnote Group

References

Untitled section

References

Associated Data

Supplementary Materials

Additional materials