Virtual screening and molecular dynamics simulations of cannabinoid receptor 2 agonists as drug candidates for osteoarthritis therapy

Abstract

Osteoarthritis is a chronic, degenerative joint disease characterized by cartilage breakdown, inflammation, and pain, significantly affecting the quality of life of affected individuals. Current treatments focus on symptom relief but do not halt disease progression, highlighting the need for novel therapeutic strategies. Cannabinoid receptor type 2 (CB2R) has emerged as a promising target for this purpose due to its role in inflammation and pain modulation, with minimal psychoactive effects compared to cannabinoid receptor type 1 (CB1R). This study performed a structure-based virtual screening to identify potential CB2R ligands with high affinity and favorable pharmacokinetic properties. Molecular docking studies led to the selection of compound P415 (PubChem CID: 154468691, 4-[[(1R,2S,6S,14R,15R,16R)-11,15-dimethoxy-5-methyl-13-oxahexacyclo [13.2.2.12,8.01,6.02,14.012,20]icosa-8(20),9,11-trien-16-yl]methoxymethyl]-3,5-dimethyl-1,2-oxazole) for further evaluation, using WIN55,212-2 (PubChem CID: 5311501, [(11R)-2-methyl-11-(morpholin-4-ylmethyl)-9-oxa-1-azatricyclo[6.3.1.0⁴,¹²]dodeca-2,4(12),5,7-tetraen-3-yl]-naphthalen-1-ylmethanone) as a reference agonist. Molecular dynamics simulations (500 ns) were also conducted in triplicate to assess the stability and dynamic behavior of CB2R-ligand complexes. The RMSD and RMSF analyses revealed distinct conformational stability patterns, with the CB2R-P415 complex showing a more stabilized profile over time with an average RMSD of 0.4 nm. Analysis of the most stable replicates through MM/PBSA binding free energy calculations yielded −51.00 kcal/mol for CB2R-P415 and -55.65 kcal/mol for CB2R-WIN55,212-2. ADMET predictions indicated that P415 possesses drug-like properties, including a log P of 3.36, TPSA of 62.95 Ų, and an LD₅₀ of 3.086 mg/kg, with no predicted AMES toxicity or hepatotoxicity. These results suggest that P415 is a promising CB2R ligand with high binding affinity and structural stability, warranting further experimental validation for potential therapeutic applications in osteoarthritis treatment.

1 Introduction

In recent years, numerous studies and clinical trials have highlighted the severity and growing impact of joint diseases, particularly osteoarthritis and rheumatoid arthritis (Kloppenburg and Berenbaum, 2020). Based on recent reports from the World Health Organization (WHO), the Global Burden of Disease (GBD) Study 2019, and updated projections from 2025, the epidemiological profile of joint diseases has become increasingly critical (Brennan-Olsen et al., 2017; Sun et al., 2025). According to the WHO, approximately 528 million people worldwide were living with osteoarthritis in 2019, representing a staggering 113% increase in prevalent cases since 1990 (Long et al., 2022). The Global Burden of Disease (GBD) Study 2019 further specifies that the global age-standardized prevalence rate (ASR) reached 6,348.25 per 100,000 in 2019, while more recent data from 2021 highlights that prevalent cases continue to rise significantly in middle- and low-income countries (Steinmetz et al., 2023). This burden is disproportionately distributed, with approximately 60% of cases occurring in women and 73% in individuals aged 55 years or older (Sun et al., 2025). The aging population, combined with an increase in risk factors such as obesity and type II diabetes, has driven the number of cases, revealing that traditional therapeutic approaches are insufficient to control symptoms and slow disease progression (Bryk and Starowicz, 2021). High body mass index (BMI) and metabolic syndrome are recognized as primary drivers, particularly in aging populations where metabolic functions decline and joint stress increases. In particular, the knee remains the most frequently affected joint, with a prevalence of 365 million cases, followed by the hand and hip (Long et al., 2022). Conventional treatments, including paracetamol, nonsteroidal anti-inflammatory drugs, corticosteroid injections, and tramadol, provide only modest pain relief. At the same time, their prolonged use can lead to significant adverse effects, such as gastrointestinal and cardiovascular complications, and even increased mortality (Philpott et al., 2017; Schmidt et al., 2018; Leopoldino et al., 2019; Zeng et al., 2019).

Once considered a simple cartilage wear-and-tear disease, osteoarthritis is now understood as a multifactorial condition in which inflammatory processes play a central role (Boutet et al., 2024). The release of pro-inflammatory mediators, such as interleukin-1β and tumor necrosis factor-α, stimulates macrophage activity and the production of degradative enzymes, promoting a self-perpetuating cycle of inflammation and joint tissue destruction (Bryk and Starowicz, 2021). This imbalance results in chronic pain, stiffness, swelling, and limited mobility, compromising patients' quality of life and generating considerable economic and social impacts. Moreover, approximately 30% of individuals with osteoarthritis exhibit characteristics of neuropathic pain, reinforcing the need for therapies that comprehensively address both inflammation reduction and pain signal modulation (Philpott et al., 2017).

Given these limitations, the search for innovative therapeutic alternatives has driven scientific efforts toward new approaches. One emerging area of interest is the study of the endocannabinoid system, a complex signaling network that plays a crucial role in regulating pain, inflammation, and immune processes (Hauser et al., 2017; Philpott et al., 2017). This system comprises endogenous cannabinoids, such as anandamide and 2-arachidonoylglycerol, as well as cannabinoid receptors (cannabinoid receptor 1 – CB1R and cannabinoid receptor 2 – CB2R), and enzymes responsible for their synthesis and degradation. CB1 receptors, predominantly located in the central nervous system, modulate neurotransmitter release and influence motor and emotional functions, while CB2 receptors, mainly present in immune cells and microglia, play a fundamental role in inflammatory responses and the maintenance of tissue integrity (Hryhorowicz et al., 2019). Recent research has shown that in degenerated joint tissues, there is an increased expression of cannabinoid receptors, particularly CB2, suggesting the existence of a compensatory mechanism aimed at reducing inflammation and promoting tissue repair (Munro et al., 1993; Sophocleous et al., 2011; Pajak et al., 2017; Shahbazi et al., 2020). This discovery has fueled studies investigating the therapeutic potential of cannabinoid modulators for the treatment of osteoarthritis (Blake et al., 2006; Fukuda et al., 2014).

Cannabis sativa, in particular its phytocannabinoids, offers potential for the development of new therapeutic interventions. While tetrahydrocannabinol (THC) is the primary psychoactive component, cannabidiol (CBD) stands out for its anti-inflammatory and immunomodulatory properties without the psychoactive effects associated with THC (Gonen and Amital, 2020). Cannabidiol acts through complex mechanisms, modulating cannabinoid receptor signaling and other systems, such as G-protein-coupled receptors and the transient receptor potential vanilloid, presenting efficacy in slowing arthritis progression and alleviating symptoms (Blake et al., 2006). In animal models, cannabidiol administration has been associated with decreased levels of inflammatory cytokines and reduced nerve fiber sensitization, contributing to pain relief and improved joint function (Malfait et al., 2000; Thomas et al., 2007; Krustev et al., 2017). These findings suggest that integrating endocannabinoid system modulators may offer an innovative therapeutic approach, providing dual benefits: controlling inflammation and promoting tissue repair. Such a strategy could represent an advantageous alternative to conventional treatments, which have limited efficacy and significant health risks for patients (Shahbazi et al., 2020).

In this context, it is essential to develop strategies that explore selective CB2R-targeting compounds, considering their structural similarity to CB1R and the need to avoid unwanted psychoactive effects (Fukuda et al., 2014). Therefore, this study aims to employ structure-based molecular modeling techniques, in particular virtual screening, to identify potential selective CB2R agonists as new therapeutic alternatives for osteoarthritis (Mizera et al., 2020; Aviz-Amador et al., 2021). By utilizing these computational tools, promising compounds were filtered and prioritized, enabling a more effective direction for future experimental assays. This work, therefore, contributes to the pursuit of innovative, safer treatments aligned with the growing need for more effective therapies with lower risks of adverse effects.

2 Materials and methods

2.1 Pharmacophore generation, validation, and structure-based filtering

A structure-based pharmacophore model was generated using the Pharmit platform (http://pharmit.csb.pitt.edu/) (Sunseri and Koes, 2016). The complete multi-stage virtual screening workflow, including the specific number of compounds evaluated at each step, is summarized in Figure 1. The co-crystallized structure of the CB2R (PDB ID: 6PT0), complexed with the agonist WIN55,212-2 (PubChem CID: 5311501, chemical name: [(11R)-2-methyl-11-(morpholin-4-ylmethyl)-9-oxa-1-azatricyclo[6.3.1.0⁴,¹²]dodeca-2,4(12),5,7-tetraen-3-yl]-naphthalen-1-ylmethanone; ligand ID: WI5), was utilized as the structural template (Xing et al., 2020). In the Pharmit interface, the PDB structure was loaded, the WI5 was selected, and the binding site was defined. Crystallographic water molecules within the binding pocket were excluded from the analysis to focus exclusively on direct receptor-ligand interactions. Subsequently, the precise pharmacophoric features were defined by uploading a manually prepared ligand structure into the “Load Features” section. This ligand file included explicitly added hydrogen atoms and atomic partial charges calculated using the semi-empirical PM7 method (MOPAC 2016), ensuring a chemically accurate representation for feature extraction (MOPAC 2016, http://OpenMOPAC.net/, Stewart Computational Chemistry, 2016). The Pharmit server was employed to analyze the geometric and chemical environment of the ligand’s pose within the binding pocket. This process identified seven features: two hydrogen bond acceptors, one aromatic region, and four hydrophobic regions, each represented by a spherical exclusion region (Figure 2). These features correspond to spatial regions where interactions are most likely to contribute to CB2 activation and were selected as the basis for subsequent virtual screening.

FIGURE 1

FIGURE 2

To ensure that the pharmacophore model was robust and could discriminate active from inactive compounds, an in silico validation step was performed using a curated dataset of 2,750 CB2R ligands from ChEMBL and BindingDB (Liu et al., 2007; Zdrazil et al., 2024). This dataset was prepared by removing duplicates, applying Pan-Assay Interference Compounds (PAINS) filters, excluding rule-of-five violators, and addressing inconsistent activity records. Ligands were classified as active (pK_i ≥ 6.5) or inactive (pK_i < 6.5)(Lipinski et al., 1997; Baell and Holloway, 2010). This curated set was uploaded to the Pharmit server as a private chemical library in SMILES format, yielding 2,708 accepted molecules and 40,523 distinct conformers. The generated pharmacophore model was applied to screen this private library, identifying active molecules based on the geometric fit score. This validation step successfully located the reference agonist, WIN55,212-2 (PubChem5311501, CHEMBL188 or BindingDB-50313648), indicating the model’s capacity to recognize known agonists (Supplementary Figure S1).

Following validation, the pharmacophore model was used to screen the PubChem database, which comprises 499,442,812 conformers from 103,302,052 molecules. The Shape filtering options within the Pharmit server were configured as follows: the Inclusive shape constraint was applied to the co-crystallized ligand with a tolerance of 1 Å, ensuring at least one heavy atom of the screened compound aligned with the pharmacophore features. Conversely, the Exclusive shape constraint was applied to the receptor with a tolerance of 1 Å, ensuring that hits with heavy-atom centers within sterically forbidden regions of the binding pocket were discarded. In addition to these shape constraints, the initial filtering for drug-likeness properties was performed based on three established pharmaceutical chemistry rules: Lipinski et al. (1997), Ghose et al. (1999), and Veber et al. (2002) (Lipinski et al., 1997; Ghose et al., 1999; Veber et al., 2002). The applied physicochemical filtering criteria were: molecular weight (MW) ≤ 500 Da, polar surface area (PSA) ≤ 140 Å2, rotatable bonds (RB) ≤ 10, hydrogen bond acceptors (HBA) ≤ 10, hydrogen bond donors (HBD) ≤ 5, and estimated Log P ranging from −0.4 to +5.6. Hits satisfying these criteria were ranked according to Pharmit’s default scoring function.

Finally, only hits meeting the following criteria were retained for further evaluation: (i) AutoDock Vina score ≤ −8 kcal/mol; (ii) minimized values of Root Mean Square Deviation (RMSD) relative to the co-crystallized ligand pose ≤2 Å; (iii) availability of a single low-energy conformer after minimization. This multi-stage filtering strategy ensured that only geometrically compatible, energetically favorable, and drug-like CB2R candidates progressed to molecular docking and molecular dynamics simulations.

2.2 Molecular docking

Molecular docking studies were conducted using GOLD 2021.1.0 software (Genetic Optimization for Ligand Docking, CCDC Software Ltd.). The docking protocol was validated through a self-docking (redocking) procedure, in which the co-crystallized agonist (WIN55,212-2) was removed and docked back into the CB2R binding site (PDB ID: 6PT0). The optimized parameters were determined by systematically evaluating: 1) different scoring functions (ASP, ChemScore, GoldScore, and ChemPLP); 2) varying binding site diameters; 3) the inclusion of structural water molecules; and 4) amino acid side-chain flexibility. The parameter combination that yielded the smallest Root Mean Square Deviation (RMSD) and the highest fitness score was selected for subsequent screening. Using this validated protocol, docking simulations were performed for the 508 molecules previously selected from the rigid docking stage in Pharmit, and, to ensure consistency and clarity throughout the manuscript, these compounds were designated P1 through P508. Atomic partial charges for these ligands were calculated using MOPAC2016 with the semi-empirical PM7 method. The ChemScore function was applied as a rescore to calculate relative binding energy data (Souza et al., 2025).

Following the docking simulations, the top 20 molecules were selected based on the lowest estimated binding energies (ΔGChemScore) and a consensus analysis of their interaction profiles relative to the reference ligand. This selection process was automated using in-house Python scripts that integrated results from the GOLD docking runs and ranked compounds based on rescoring performance. The rationale for selecting the 20 best candidates balanced computational efficiency with a sufficiently diverse chemical space for subsequent ADMET analyses and molecular dynamics simulations.

2.3 In silico ADMET (absorption, distribution, metabolism, excretion, and toxicity) prediction

The design of bioactive molecules also requires assessing their ADME properties to understand their physicochemical behavior in the human body. Moreover, the toxicity profile (T) of these molecules and their metabolites is also essential for evaluating their potential adverse effects on biological systems (Daoui et al., 2022). The Lipinski’s Rule of Five was applied using the SwissADME server (Daina et al., 2017). The SMILES code for the 20 selected molecules was used to predict their physicochemical properties and generate a Boiled-Egg graph, which provides insight into gastrointestinal absorption (GA) and blood-brain barrier (BBB) permeation (Daina and Zoete, 2016). For the subsequent analyses, molecules predicted to cross the blood-brain barrier (BBB) were prioritized. This selection criterion is based on the functional expression of CB2R within the central nervous system (CNS) and its role in modulating neuroinflammation and central sensitization. Given the neuropathic component of osteoarthritis pain, identifying ligands capable of modulating both the central and peripheral immune systems was deemed essential for achieving comprehensive therapeutic efficacy (Hryhorowicz et al., 2019). Additionally, toxicity prediction was performed using the pkCSM server (https://biosig.lab.uq.edu.au/pkcsm/), yielding results for the Ames test, hepatotoxicity, and skin permeability (Pires et al., 2015).

2.4 Selectivity test by ALPACA

The promising molecules were submitted for further affinity and selectivity analysis using the ALPACA platform (“A machine Learning Platform for Affinity and selectivity profiling of CAnnabinoid receptors modulators”), which enables the use of machine learning-based affinity classifiers for CB1 and CB2 receptors, as well as CB2/CB1 selectivity (https://www.ba.ic.cnr.it/softwareic/alpaca/) (Delre et al., 2023). Thus, in addition to predicting potential interactions and bioavailability, molecules with predicted affinity and selectivity were selected.

2.5 Molecular dynamics simulations

The selected candidate, P415 (PubChem CID: 154468691; chemical name: 4-[[(1R,2S,6S,14R,15R,16R)-11,15-dimethoxy-5-methyl-13-oxahexacyclo [13.2.2.12,8.01,6.02,14.012,20]icosa-8(20),9,11-trien-16-yl]methoxymethyl]-3,5-dimethyl-1,2-oxazole), and the reference agonist, WIN55,212-2, were subjected to molecular dynamics simulations. The structures of the ligand-receptor complexes in their membrane environments were mapped to coarse-grained (CG) models. To add POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) bilayer lipids, MemProtMD was used, while the membrane-spanning regions and protein orientation were predicted using Memembed (Stansfeld et al., 2015). The protein structure was converted into a CG representation, embedded in a lipid membrane, and hydrated with water and ions to balance the system using martinize.py v2.5 and Insane from the Martini repository. Molecular dynamics simulations were performed for 10 ns using GROMACS, and the final snapshot was back-mapped to the CHARMM36m representation using CG2AT (Huang et al., 2016; Vickery and Stansfeld, 2021; Bauer et al., 2022).

After the system preparation, the all-atom simulation was conducted in three stages, using a force constant of 1,000 kJ·mol^-1·nm^-2 to constrain all heavy atoms. The first stage, energy minimization, involved an initial optimization of each system geometry using the steepest-descent algorithm over 5,000 iterations (5 ps). The subsequent step involved system equilibration in two stages, each running for 100,000 iterations (10 ns). The first equilibration stage was performed at constant particle number, volume, and temperature (NVT ensemble). The second stage was carried out at constant particle number, pressure, and temperature (NPT ensemble) at 1 atm and 310K, employing the C-rescale barostat and the V-rescale thermostat, respectively (Bernetti and Bussi, 2020). The Particle Mesh Ewald (PME) algorithm was used to calculate electrostatic interactions, and bonds were constrained using the LINCS algorithm (Hess et al., 1997).

Finally, a 500 ns production run under constant pressure was performed, in triplicate. Comparative data calculations, including RMSD and Root Mean Square Fluctuation (RMSF), were performed using GROMACS’s integrated tools to analyze MD trajectories.

Ultimately, the binding free energy between ligands and proteins was estimated using the GROMACS module “gmx_MMPBSA” (Kumari et al., 2014; Valdés-Tresanco et al., 2021). Additionally, an interaction map was generated using in-house Python scripts that employed the pandas, matplotlib, and seaborn libraries, along with the Protein-Ligand Interaction Profiler (PLIP) 2.4.0 (Salentin et al., 2015). PLIP was used to extract molecular interaction data from the simulations (De Paula et al., 2025). These data were compiled into a single dataframe using Pandas and visualized as scatter plots, with distinct colors for different interaction types using Matplotlib and Seaborn.

3 Results

3.1 Pharmacophore model generation and preliminary evaluation

The study was initiated by constructing a structure-based pharmacophore model using known CB2 receptor agonists. This model served as the foundation for the virtual screening of compounds from the PubChem database. Table 1 presents the features of the generated pharmacophore, detailing the types and number of identified pharmacophoric features.

To identify potential agonists of the CB2R (PDB ID: 6PT0), 499,442,812 conformers were screened from 103,302,052 molecules in the PubChem database using a pharmacophore model (Table 1). This process resulted in 1,560 small-molecule hits, which were further filtered by rigid docking using the Vina interface on the Pharmit server, applying an RMSD cutoff of 2 Å and a maximum Vina score of ≤ −8 kcal/mol. This process yielded 508 small-molecule hits for further docking simulations using the GOLD program. The redocking validation step confirmed the robustness of the protocol, returning the best configurations as a function of the ASP score, a 10 Å site diameter, and ChemScore rescoring. This combination resulted in an RMSD value of 0.8494 Å and a Biggest Fitness Score (BFS) of 54.6357 (Supplementary Table S1). The low RMSD value (<1.0 Å) indicates that the docking protocol was highly accurate in replicating the experimental binding pose (Supplementary Table S1). Afterwards, 20 compounds were selected as the most promising candidates (Table 2) since they had the lowest predicted binding energies (ΔG_ChemScore).

TABLE 2

PubChem CID	ID	Interactions	*ΔG (kcal/mol)	PubChem CID	ID	Interactions	*ΔG
155092381	P417		−53.18	6825440	P40		−54.47
1252403	P31		−53.22	89715247	P01		−54.61
23804736	P64		−53.27	10093112	P22		−54.61
118194531	P02		−53.57	15838126	P28		−54.78
8094022	P46		−53.87	159612009	P29		−54.82
4527323	P66		−53.9	154468691	P415		−54.87
10454310	P136		−54.14	19588151	P101		−54.93
68033057	P08		−54.23	88094023	P44		−55.91
117807624	P88		−54.27	145609961	P25		−55.99
11743324	P35		−54.45	124008102	P114		−56.54
5311501	WIN55,212-2		−53.07	​

20 compounds with the lowest binding energy (ΔG_ChemScore) values as the most promising candidates and 2D representation of interactions.

*

ΔG: values calculated by the ChemScore function (kcal/mol). Interactions: hydrogen bonding (dark green), van der Waals (light green), π stacking (pink), saline bridge (orange), halogen (cyan).

It is important to mention that all the compounds share at least two types of similar interactions with the crystallographic ligand, such as π - π interactions with the catalytic residues (Phe87, Phe90, Phe94, Ile110, Val113, Phe117, Phe183, Pro184, Trp194, Trp258, Val261, and Phe281), as well as hydrogen bonding with His95 and Ser285. These interactions suggest that these compounds may bind to the CB2R site with stability, but further validation by molecular dynamics simulations is needed.

3.2 ADMET properties

Predicting pharmacokinetic and toxicological properties is essential for prioritizing promising compounds in the early stages of drug development. Analyses using SwissADME, pKCSM, and ALPACA enabled assessment of the selected compounds' permeability, bioavailability, and potential toxicological risks. The pharmacokinetic analysis revealed that the studied compounds exhibit distinct profiles in terms of absorption and permeability. The compounds' projections in the WLOGP vs. TPSA (Boiled-Egg) plot (Figure 3) were used to predict their absorption in the human body.

FIGURE 3

Figure 3 shows three regions: yellow (egg yolk), white (egg white), and gray. The yellow region indicates compounds with a high likelihood of crossing the blood-brain barrier; the white region indicates compounds likely to be absorbed by the gastrointestinal tract; and the gray region indicates compounds unlikely to be absorbed in any form. This plot indicates that most candidates are located in regions compatible with permeability across the BBB, suggesting potential action in the central nervous system. However, compounds outside this region may have a lower capacity to cross the BBB, limiting their efficacy. In this analysis, a few molecules (including WIN55,212-2) were found in the yellow region. However, only P66 (PubChem CID: 4527323, chemical name: N-[9,10-dioxo-4-(2,4,6-trimethylanilino)anthracen-1-yl]-N-methylacetamide), P64 (PubChem CID: 23804736, chemical name: 14-methyl-10-(2,4,6-trimethylanilino)-14-azatetracyclo[7.7.1.02,7.013,17]heptadeca-1(16),2,4,6,9(17),10-hexaene-8,15-dione), and P415 were predicted, including PGP (non-substrate of P-gp). This result supports our goal of identifying a selective compound that is less likely to be actively pumped out of cells by this efflux transporter. Therefore, P415 was selected for further analysis by molecular dynamics simulations.

In addition, the potential toxicity of the selected compounds was also evaluated, including AMES toxicity and hepatotoxicity predictions. These factors should be considered when screening the most promising candidates for further experimental studies. Compounds P64 and P66 presented a more favorable profile, with no predicted hepatotoxicity. Despite the hepatotoxicity prediction for P415, its predicted LD₅₀ is the highest one, suggesting that the proposed administered dose could be flexible.

Selectivity for the CB2 receptor over CB1R is a crucial criterion for minimizing adverse effects associated with CB1R activation, such as sedation and psychotropic effects. The predicted pK_i values (Table 3) indicate that compounds P64, P66, and P415 have significant affinity for CB2R, with values above 6.5, making them viable candidates for further investigation. However, only P415 maintains a pK_i above 7, suggesting that it interacts best with the target.

3.3 Molecular dynamics simulations

To evaluate the stability of the complexes formed by CB2R-WIN55,212-2 and CB2R-P415 over time, molecular dynamics simulations (500 ns) were performed. Initially, the variation in the RMSD and RMSF values of residues were analysed. RMSD measures the deviation of the protein structure from its initial to its final conformation. The RMSD values obtained during the simulations reflect the protein’s structural stability. The RMSD plots of the ligand-protein complexes showed low RMSD values, indicating that the systems were well-equilibrated and stable over time (Figure 4).

FIGURE 4

The RMSD plot (Figure 4A) shows that the apo CB2 receptor exhibited the least conformational variation over time, with values near 0.3 nm. The CB2R-WIN55,212-2 complex showed a progressive increase in RMSD and greater variation throughout the simulation, stabilizing at approximately 0.4 nm after 300 ns. Meanwhile, the CB2R-P415 complex exhibited an initial variation of the RMSD values until 100 ns, reaching an average value of 0.4 nm, and then remained stable until the end of the simulation, suggesting a more significant impact of ligand interaction on the receptor stability (Supplementary Figure S2).

The RMSF analysis (Figure 4B) reveals that structural fluctuations are similar among the systems, with a prominent peak in the range related to residues 227-239, indicating a more flexible region associated with ICL3 (intracellular loop 3), which is structurally linked to receptor interaction with the G protein (structural regions were identified in Supplementary Figure S3) (Shao et al., 2016; Li et al., 2019). The CB2R-P415 complex exhibits slightly lower fluctuations than the apo form and the reference agonist complex, with only the ECL1 region (extracellular loop 1) showing slightly more pronounced fluctuations. This may indicate a differential impact on the receptor’s local dynamics compared to CB2R-WIN55,212-2. These results suggest that although both compounds interact with CB2R, the receptor’s conformational stability may be modulated differently, with the CB2R-P415 complex potentially stabilizing over time.

An interaction analysis was also performed, comparing the duration of interactions during the molecular dynamics (MD) simulations of the complexes with P415 and WIN55,212-2 (Figure 5). Over time, in the CB2R-P415 complex, the ligand interacted with at least 7 residues, similar to those observed in interactions with CB2R-WIN55-212,2. These residues are Tyr25, Ser90, His95, Thr114, Phe117, Phe281 and Ser285. Thus, the interaction patterns in both complexes were compared. A combination of hydrophobic and π-π interactions stabilizes the CB2R-P415 and CB2R-WIN55-212,2 complexes. In the CB2R-P415 complex, stability is predominantly conferred by hydrogen bonds, mainly with Ile27, Thr114, Ala282, and Ser285 in addition to π-stacking interactions with Phe87, Phe91, and Phe117. In comparison, the reference complex CB2R-WIN55,212-2 is stabilized by π-stacking interactions with Tyr25, Phe87, Phe94, Phe117, and Phe281, complemented by hydrogen bonds with Ser90, Thr114, and Ser285. A critical distinction lies in the π-cation interaction involving the Phe183 residue and the WIN55,212-2 ligand, which provides superior binding strength compared to non-cationic interactions due to strong electrostatic attraction.

FIGURE 5

The interaction affinity is usually defined as the change in free energy upon formation of a ligand-receptor complex. This parameter is directly linked to the ligand’s binding strength. The results show that the CB2R-P415 complex has the lowest binding free energy, indicating a more energetically favorable interaction (Table 4). This value is mainly derived from contributions from π-π and hydrogen-bond interactions, which were also observed in other studies, in which stability at the catalytic site is generally associated with contributions from ΔE_vdw (Liu and Liu, 2019; Wang et al., 2020). Due to the predominance of aromatic rings in the selected compound, the electronic delocalization near the aromatic catalytic residues also favored this type of contribution, potentially strengthening the binding of P415 and stabilizing the ligand in the binding site.

Moreover, when analyzing the energy graph over time (Figure 6A), the energy fluctuations for the CB2R-P415 complex exhibit oscillations with an average energy value around −51.00 kcal/mol, but with peaks reaching more negative values, close to −60 kcal/mol. These fluctuations are characteristic of the dynamic interactions between the receptor and the ligand during the simulation. The complex’s stability appears reasonable. However, some temporal variations may reflect the receptor’s adaptability or flexibility during ligand binding. These results indicate the overall stabilization of the CB2R-P415 complex and suggest a potential agonist effect.

FIGURE 6

In the case of the CB2R-WIN55-212,2 complex, it exhibits similar behavior but with more constant energy throughout the simulation. Fluctuations are also present, but they seem less pronounced than those observed in the CB2R-P415 complex, suggesting lower instability or greater stability in the interaction. Further detailed analysis of the specific energy contributions of different residues in P415 and WIN55-212,2 can be conducted from the plots in Figures 6B,C (more details of the energy decomposition per residue are displayed in Supplementary Figure S3). In the CB2R-P415 complex, the residue Phe87 made a significant contribution, possibly playing an anchoring role in stabilizing the ligand at the binding site. In the CB2R-WIN55-212,2 complex, the π-cation interaction involving Phe183 and the WIN55-212,2 ligand is a distinct type of interaction, which has a critical importance in the stability of the complex and the binding affinity, given its strength exceeding that of π-stacking interactions. This is due to the attraction between the aromatic ring and the positive charge on the ligand, which provides a more stable bond, as the electrostatic interaction between the π orbitals of the aromatic ring and the cation in the ligand is stronger than that between two aromatic rings.

4 Discussion

This study employed a robust, multi-faceted approach to identify novel, potent, and selective CB2R agonists, a key receptor in the endocannabinoid system (ECS). Given the pivotal role of CB2R in various physiological processes, including immune response modulation and inflammation, targeting this receptor holds significant therapeutic promise (Alfei et al., 2023). The computational workflow pipeline employed in this study, integrating pharmacophore modeling, virtual screening, molecular docking, in silico ADMET prediction, molecular dynamics simulations, and MM-PBSA end-state free energy calculations, is well-supported by recent literature focusing on the discovery of novel CB2R ligands, to comprehensively evaluate the possible efficacy, selectivity, and binding affinity of candidate compounds (Uba et al., 2022). Using the crystal structure of WIN55,212-2 bound to CB2R as a reference, novel compounds with improved potency and selectivity for CB2R were designed and identified (Mohammadi Vosough et al., 2019; Xing et al., 2020).

In this study, 20 compounds were identified by molecular docking and exhibited favorable interactions with key catalytic residues of CB2R, as indicated by estimated ΔG values. These compounds presented a combination of hydrophobic, hydrogen bonds, and π-stacking interactions with critical residues in the receptor’s binding pocket. Among these compounds, P415 stood out for its promising physicochemical properties and alignment with Lipinski’s rule of five, a criterion for drug-like properties. The Boiled-Egg plot indicated that P415 has the potential to cross the blood-brain barrier, an essential characteristic for therapeutic efficacy in neurological disorders. This feature is significant, given the increasing interest in ECS as a therapeutic target for brain-related diseases (Hauser et al., 2017; Uba et al., 2022).

The toxicity profile of P415 also supports its further development as a lead compound. In silico toxicity predictions showed no significant concerns, which is crucial for advancing the compound to in vitro and in vivo testing stages (Ekins et al., 2000). Moreover, the molecular dynamics simulations provided significant insights into the stability and dynamic behavior of the CB2R-P415 complex. The low RMSD values observed in these simulations suggest that the complex remains stable throughout the simulation, reflecting a well-defined and stable binding interaction between the ligand and receptor (Aviz-Amador et al., 2021). This stability is crucial for maintaining the receptor’s function and preventing ligand dissociation, which could compromise therapeutic efficacy (Bahrami et al., 2018).

The RMSF analysis further revealed the receptor’s dynamic nature, particularly at specific residues involved in ligand binding. Distinct mobility profiles in the intra- and extracellular loop regions were observed, with the highest fluctuations occurring in the extracellular loop region, a behavior that corroborates the findings of Feng et al. (2014). In their dynamic model of CB2R, the authors highlighted that ECL2 and ICL3 exhibit significant intrinsic plasticity, which is essential for accommodating ligand entry and facilitating G-protein coupling, respectively. In the CB2R-P415 complex, the magnitude of these fluctuations was comparable to those reported in the literature for agonists, suggesting that P415 binding maintains the receptor’s structural integrity while allowing the flexibility required for activation (Feng et al., 2014). In particular, the relative stabilization observed in ECL2 during the 500 ns simulation indicates a robust interaction network anchoring the ligand, resembling the dynamic stability patterns. Fluctuations observed suggest that the receptor undergoes conformational changes upon ligand binding, which may help stabilize the complex. These findings align with those of Rao et al. (2023), who highlighted the importance of protein flexibility in receptor-ligand interactions. The ability of P415 to stabilize the receptor while also allowing for the necessary conformational changes suggests its potential as a stable and effective agonist (Hua et al., 2020).

Furthermore, the stability of the CB2R-P415 complex was evaluated considering the conformational state of the “toggle switch” residue, Trp258. This residue is a critical molecular sensor in Class A GPCRs, where its side-chain orientation directly correlates with the transition between inactive and active receptor states. In alignment with other findings, the maintenance of specific interactions, the stabilization of Trp258 throughout molecular docking, and the 500 ns trajectory suggest that P415 acts as a potential agonist, capable of stabilizing the active-like conformation required for signal transduction without the off-target effects associated with CB1R activation (Uba et al., 2022; Bousselham et al., 2025). The predicted free binding energy further confirmed the strength of the interaction between P415 and CB2R. Negative values of the free binding energy indicate a highly favorable binding interaction with P415, suggesting a higher binding affinity than that of other candidates. This is primarily due to specific interactions with residues such as Phe87 and Ala282, which contribute significantly to the complex’s stability. Interestingly, when comparing P415 to other ligands, such as those interacting through π-stacking (like Phe183 in the CB2R-WIN55-212,2 complex), the presence of π-cation interactions with WIN55-212,2 seems to provide a more robust and stable complex, making WIN55-212,2 potentially more stable in its binding than other non-cationic ligands (Omar et al., 2022; Bousselham et al., 2025). However, the overall negative free binding energy related to P415 still suggests that it is a potent and stable ligand for CB2R. The use of MM-PBSA as a rescoring tool provides a more accurate estimate of binding affinity than standard docking scores, accounting for receptor flexibility and solvent effects, ensuring that the final lead candidate, P415, possesses a robust binding profile (Uba et al., 2022).

These findings are significant as they highlight P415 as a promising candidate for drug development targeting CB2R. Its unique properties could offer an advantage in therapeutic strategies, particularly in diseases related to endocannabinoid system (ECS) dysregulation (Krustev et al., 2017). The stability and affinity of P415, with its favorable physicochemical properties and safety profile, make it an ideal candidate for further preclinical testing. Ultimately, these results underscore the potential of combining various computational techniques, such as virtual screening, molecular docking, and molecular dynamics simulations, to inform subsequent experimental studies and expedite the discovery of novel drug candidates targeting ECS.

5 Limitations and future prospects

Although the computational strategies employed in this study (including molecular docking, molecular dynamics, and MM-PBSA calculations) provided robust evidence for P415’s potential as a selective CB2R agonist, certain limitations must be acknowledged. Firstly, the findings are strictly based on in silico models, which, despite their high predictive power and grounding in established crystallographic structures, require experimental validation. Subsequent in vitro assays, such as radioligand binding studies and cAMP functional assays, are necessary to confirm the binding affinity and the specific agonist efficacy of the identified compounds. Furthermore, while the ADMET profile suggests favorable pharmacokinetics and low toxicity, in vivo studies will be essential to evaluate the actual permeability of the blood-brain barrier and the long-term safety in complex biological systems. Prospects include structural optimization of the P415 scaffold to further enhance its selectivity for CB1R and the exploration of its therapeutic impact in specific animal models of osteoarthritis to bridge the gap between computational discovery and clinical application.

6 Conclusion

In this study, a multi-stage virtual screening workflow was successfully implemented to identify novel selective CB2R agonists with potential applications in the treatment of osteoarthritis. From the integration of pharmacophore modeling, molecular docking, and 500 ns molecular dynamics simulations, compound P415 was identified as a promising lead. The results indicated that P415 maintains a stable interaction network within the CB2R binding pocket, characterized by favorable binding free energy values and consistent stabilization of key structural regions, including the intra- and extracellular loops. The energetic superiority of the reference ligand WIN55,212-2, driven by specific π-cation interactions, provided a benchmark for understanding the binding efficiency of the newly identified scaffold. Moreover, the predicted ADMET profile and the high LD₅₀ values suggest that P415 is a safe candidate for further development, potentially offering a therapeutic alternative that addresses both inflammatory and neuropathic pain without the psychoactive risks associated with CB1R activation. Therefore, this study contributes to the growing field of endocannabinoid modulation and provides a solid foundation for future experimental investigations into CB2R-targeted therapies.

References

• 1

  AlfeiS.SchitoG. C.SchitoA. M. (2023). Synthetic pathways to non-psychotropic phytocannabinoids as promising molecules to develop novel antibiotics: a review. Pharmaceutics15 (7), 1889. 10.3390/PHARMACEUTICS15071889

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  Aviz-AmadorA.Contreras-PuentesN.Mercado-CamargoJ. (2021). Virtual screening using docking and molecular dynamics of cannabinoid analogs against CB1 and CB2 receptors. Comput. Biol. Chem.95, 107590. 10.1016/J.COMPBIOLCHEM.2021.107590

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BaellJ. B.HollowayG. A. (2010). New substructure filters for removal of pan assay interference compounds (PAINS) from screening libraries and for their exclusion in bioassays. J. Medicinal Chemistry53 (7), 2719–2740. 10.1021/jm901137j

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BahramiH.SalehabadiH.NazariZ.AmanlouM. (2018). Combined virtual screening, DFT calculations and molecular dynamics simulations to discovery of potent MMP-9 inhibitors. Lett. Drug Des. and Discov.16 (8), 892–903. 10.2174/1570180815666181008095950

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BauerP.HessB.LindahlE. (2022). GROMACS 2022.3 source code. 10.5281/ZENODO.7037338

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BernettiM.BussiG. (2020). Pressure control using stochastic cell rescaling. J. Chem. Phys.153 (11), 114107. 10.1063/5.0020514/199610

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BlakeD. R.RobsonP.HoM.JubbR. W.McCabeC. S. (2006). Preliminary assessment of the efficacy, tolerability and safety of a cannabis-based medicine (sativex) in the treatment of pain caused by rheumatoid arthritis. Rheumatology45 (1), 50–52. 10.1093/RHEUMATOLOGY/KEI183

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BousselhamF.MouhcineM.GhichaI.KadilY.BaghrousS.MoundA.et al (2025). Cannabis compounds: docking and dynamics study. J. Drug Deliv. and Ther.15 (9), 83–91. 10.22270/jddt.v15i9.7371

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BoutetM. A.NervianiA.Fossati-JimackL.Hands-GreenwoodR.AhmedM.RivelleseF.et al (2024). Comparative analysis of late-stage rheumatoid arthritis and osteoarthritis reveals shared histopathological features. Osteoarthr. Cartil.32 (2), 166–176. 10.1016/j.joca.2023.10.009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  Brennan-OlsenS. L.CookS.LeechM. T.BoweS. J.KowalP.NaidooN.et al (2017). Prevalence of arthritis according to age, sex and socioeconomic status in six low and middle income countries: analysis of data from the world health organization study on global AGEing and adult health (SAGE) wave 1. BMC Musculoskeletal Disorders18 (1), 271. 10.1186/s12891-017-1624-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  BrykM.StarowiczK. (2021). Cannabinoid-based therapy as a future for joint degeneration. Focus on the role of CB2 receptor in the arthritis progression and pain: an updated review. Pharmacol. Rep.73 (3), 681–699. 10.1007/S43440-021-00270-Y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  DainaA.ZoeteV. (2016). A BOILED‐Egg to predict gastrointestinal absorption and brain penetration of small molecules. Chemmedchem11 (11), 1117. 10.1002/CMDC.201600182

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  DainaA.MichielinO.ZoeteV. (2017). SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Reports7 (1), 1–13. 10.1038/srep42717

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  DaouiO.MazoirN.BakhouchM.SalahM.BenharrefA.Gonzalez-ColomaA.et al (2022). 3D-QSAR, ADME-Tox, and molecular docking of semisynthetic triterpene derivatives as antibacterial and insecticide agents. Struct. Chem.33 (4), 1063–1084. 10.1007/s11224-022-01912-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  de PaulaH.SouzaF.FerreiraL.SilvaJ. A. B.RibeiroR.VilachãJ.et al (2025). Semisynthetic flavonoids as GSK-3β inhibitors: computational methods and enzymatic assay. Targets3 (2), 13. 10.3390/targets3020013

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  DelreP.ContinoM.AlbergaD.SavianoM.CorrieroN.MangiatordiG. F. (2023). ALPACA: a machine learning platform for affinity and selectivity profiling of CAnnabinoids receptors modulators. Comput. Biol. Med.164, 107314. 10.1016/J.COMPBIOMED.2023.107314

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  EkinsS.WallerC. L.SwaanP. W.CrucianiG.WrightonS. A.WikelJ. H. (2000). Progress in predicting human ADME parameters in silico. J. Pharmacological Toxicological Methods44 (1), 251–272. 10.1016/s1056-8719(00)00109-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  FengZ.AlqarniM. H.YangP.TongQ.ChowdhuryA.WangL.et al (2014). Modeling, molecular dynamics simulation, and mutation validation for structure of cannabinoid receptor 2 based on known crystal structures of GPCRs. J. Chemical Information Modeling54 (9), 2483–2499. 10.1021/ci5002718

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  FukudaS.KohsakaH.TakayasuA.YokoyamaW.MiyabeC.MiyabeY.et al (2014). Cannabinoid receptor 2 as a potential therapeutic target in rheumatoid arthritis. BMC Musculoskelet. Disord.15 (1), 1–10. 10.1186/1471-2474-15-275/FIGURES/5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  GhoseA. K.ViswanadhanV. N.WendoloskiJ. J. (1999). A knowledge-based approach in designing combinatorial or medicinal chemistry libraries for drug discovery. 1. A qualitative and quantitative characterization of known drug databases. J. Combinatorial Chemistry1 (1), 55–68. 10.1021/cc9800071

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  GonenT.AmitalH. (2020). Cannabis and cannabinoids in the treatment of rheumatic diseases. Rambam Maimonides Med. J.11 (1), e0007. 10.5041/RMMJ.10389

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  HauserA. S.AttwoodM. M.Rask-AndersenM.SchiöthH. B.GloriamD. E. (2017). Trends in GPCR drug discovery: new agents, targets and indications. Nat. Rev. Drug Discov.16 (12), 829–842. 10.1038/nrd.2017.178

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  HessB.BekkerH.BerendsenH. J. C.FraaijeJ. G. E. M. (1997). LINCS: a linear constraint solver for molecular simulations. J. Computational Chemistry18 (12), 1463–1472. 10.1002/(sici)1096-987x(199709)18:12<1463::aid-jcc4>3.0.co;2-h

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  HryhorowiczS.Kaczmarek-RyśM.AndrzejewskaA.StaszakK.HryhorowiczM.KorczA.et al (2019). Allosteric modulation of cannabinoid receptor 1—Current challenges and future opportunities. Int. J. Mol. Sci.20 (23), 5874. 10.3390/IJMS20235874

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  HuaT.LiX.WuL.Iliopoulos-TsoutsouvasC.WangY.WuM.et al (2020). Activation and signaling mechanism revealed by cannabinoid Receptor-Gi complex structures. Cell180 (4), 655–665.e618. 10.1016/J.CELL.2020.01.008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  HuangJ.RauscherS.NawrockiG.RanT.FeigM.De GrootB. L.et al (2016). CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat. Methods14 (1), 71–73. 10.1038/nmeth.4067

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  KloppenburgM.BerenbaumF. (2020). Osteoarthritis year in review 2019: epidemiology and therapy. Osteoarthr. Cartil.28 (3), 242–248. 10.1016/j.joca.2020.01.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  KrustevE.MuleyM. M.McDougallJ. J. (2017). Endocannabinoids inhibit neurogenic inflammation in murine joints by a non-canonical cannabinoid receptor mechanism. Neuropeptides64, 131–135. 10.1016/J.NPEP.2016.08.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  KumariR.KumarR.LynnA. (2014). G-mmpbsa -A GROMACS tool for high-throughput MM-PBSA calculations. J. Chem. Inf. Model.54 (7), 1951–1962. 10.1021/CI500020M/SUPPL_FILE/CI500020M_SI_001.PDF

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  LeopoldinoA. O.MacHadoG. C.FerreiraP. H.PinheiroM. B.DayR.McLachlanA. J.et al (2019). Paracetamol versus placebo for knee and hip osteoarthritis. Cochrane Database Systematic Reviews2 (2), CD013273. 10.1002/14651858.CD013273

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  LiX.HuaT.VemuriK.HoJ. H.WuY.WuL.et al (2019). Crystal structure of the human cannabinoid receptor CB2. Cell176 (3), 459–467.e413. 10.1016/j.cell.2018.12.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  LipinskiC. A.LombardoF.DominyB. W.FeeneyP. J. (1997). Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv. Drug Deliv. Rev.23 (1-3), 3–25. 10.1016/S0169-409X(96)00423-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  LiuY.LiuS.-Y. (2019). Exploring the strength of a hydrogen bond as a function of steric environment using 1, 2-azaborine ligands and engineered T4 lysozyme receptors. Org. and Biomolecular Chemistry17 (29), 7002–7006. 10.1039/c9ob01008d

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  LiuT.LinY.WenX.JorissenR. N.GilsonM. K. (2007). BindingDB: a web-accessible database of experimentally determined protein–ligand binding affinities. Nucleic Acids Research35 (Suppl. l_1), D198–D201. 10.1093/nar/gkl999

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  LongH.LiuQ.YinH.WangK.DiaoN.ZhangY.et al (2022). Prevalence trends of site‐specific osteoarthritis from 1990 to 2019: findings from the global burden of disease study 2019. Arthritis and Rheumatology74 (7), 1172–1183. 10.1002/art.42089

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  MalfaitA. M.GallilyR.SumariwallaP. F.MalikA. S.AndreakosE.MechoulamR.et al (2000). The nonpsychoactive cannabis constituent cannabidiol is an oral anti-arthritic therapeutic in murine collagen-induced arthritis. Proc. Natl. Acad. Sci. U. S. A.97 (17), 9561–9566. 10.1073/PNAS.160105897

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  MizeraM.LatekD.Cielecka-PiontekJ. (2020). Virtual screening of C. sativa constituents for the identification of selective ligands for cannabinoid receptor 2. Int. J. Mol. Sci.21 (15), 5308. 10.3390/IJMS21155308

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  Mohammadi VosoughE.Baradaran RahimiV.MasoudS. A.MirkarimiH. R.DemnehM. K.AbedA.et al (2019). Evaluation of protective effects of non-selective cannabinoid receptor agonist WIN 55,212-2 against the nitroglycerine-induced acute and chronic animal models of migraine: a mechanistic study. Life Sci.232, 116670. 10.1016/J.LFS.2019.116670

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  MunroS.ThomasK. L.Abu-ShaarM. (1993). Molecular characterization of a peripheral receptor for cannabinoids. Nature365 (6441), 61–65. 10.1038/365061a0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  OmarA. M.AljahdaliA. S.SafoM. K.MohamedG. A.IbrahimS. R. M. (2022). Docking and molecular dynamic investigations of phenylspirodrimanes as cannabinoid receptor-2 agonists. Molecules28 (1), 44. 10.3390/molecules28010044

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  PajakA.KostrzewaM.MalekN.KorostynskiM.StarowiczK. (2017). Expression of matrix metalloproteinases and components of the endocannabinoid system in the knee joint are associated with biphasic pain progression in a rat model of osteoarthritis. J. Pain Res.10, 1973–1989. 10.2147/JPR.S132682

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  PhilpottH. T.O'BrienM.McDougallJ. J. (2017). Attenuation of early phase inflammation by cannabidiol prevents pain and nerve damage in rat osteoarthritis. Pain158 (12), 2442–2451. 10.1097/J.PAIN.0000000000001052

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  PiresD. E. V.BlundellT. L.AscherD. B. (2015). pkCSM: predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures. J. Med. Chem.58 (9), 4066–4072. 10.1021/ACS.JMEDCHEM.5B00104/SUPPL_FILE/JM5B00104_SI_001.PDF

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  RaoK. Y.BashaS. J.MonikaK.SreelakshmiM.SivakumarI.MallikarjunaG.et al (2023). Synthesis and anti-Alzheimer potential of novel \x{03b1}-amino phosphonate derivatives and probing their molecular interaction mechanism with acetylcholinesterase. Eur. J. Med. Chem.253, 115288. 10.1016/j.ejmech.2023.115288

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  SalentinS.SchreiberS.HauptV. J.AdasmeM. F.SchroederM. (2015). PLIP: fully automated protein–ligand interaction profiler. Nucleic Acids Res.43 (Web Server issue), W443. 10.1093/NAR/GKV315

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  SchmidtM.SørensenH. T.PedersenL. (2018). Diclofenac use and cardiovascular risks: series of nationwide cohort studies. BMJ Clin. Research ed.362, k3426. 10.1136/BMJ.K3426

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  ShahbaziF.GrandiV.BanerjeeA.TrantJ. F. (2020). Cannabinoids and cannabinoid receptors: the story so far. iScience23 (7), 101301. 10.1016/J.ISCI.2020.101301

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  ShaoZ.YinJ.ChapmanK.GrzemskaM.ClarkL.WangJ.et al (2016). High-resolution crystal structure of the human CB1 cannabinoid receptor. Nature540, 602–606. 10.1038/nature20613

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  SophocleousA.Landao-BassongaE.Van't HofR. J.IdrisA. I.RalstonS. H. (2011). The type 2 cannabinoid receptor regulates bone mass and ovariectomy-induced bone loss by affecting osteoblast differentiation and bone formation. Endocrinology152 (6), 2141–2149. 10.1210/EN.2010-0930

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  SouzaF. F.VilachãJ. F.CamposO. S.de PaulaH. (2025). Prediction of novel insecticides for malaria prevention: virtual screening and molecular dynamics of ag ache inhibitors. Drugs Drug Candidates4 (3), 41. 10.3390/ddc4030041

  • CrossRef
  • Google Scholar

• 51

  StansfeldP. J.GooseJ. E.CaffreyM.CarpenterE. P.ParkerJ. L.NewsteadS.et al (2015). MemProtMD: automated insertion of membrane protein structures into explicit lipid membranes. Struct. Engl.23 (7), 1350. 10.1016/J.STR.2015.05.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  SteinmetzJ. D.CulbrethG. T.HaileL. M.RaffertyQ.LoJ.FukutakiK. G.et al (2023). Global, regional, and national burden of osteoarthritis, 1990–2020 and projections to 2050: a systematic analysis for the global burden of disease study 2021. Lancet Rheumatology5 (9), e508–e522. 10.1016/S2665-9913(23)00163-7

  • CrossRef
  • Google Scholar

• 53

  SunJ.HanS.LiangJ.LiuW.XingZ.LiQ.et al (2025). The epidemiology and burden of global osteoarthritis: a systematic analysis based on the global burden of disease database. Aging Adv.2 (4), 123–131. 10.4103/agingadv.agingadv-d-25-00003

  • CrossRef
  • Google Scholar

• 54

  SunseriJ.KoesD. R. (2016). Pharmit: interactive exploration of chemical space. Nucleic Acids Res.44 (W1), W442–W448. 10.1093/NAR/GKW287

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  ThomasA.BaillieG. L.PhillipsA. M.RazdanR. K.RossR. A.PertweeR. G. (2007). Cannabidiol displays unexpectedly high potency as an antagonist of CB1 and CB2 receptor agonists in vitro. Br. Journal Pharmacology150 (5), 613–623. 10.1038/SJ.BJP.0707133

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  UbaA. I.AluwalaH.LiuH.WuC. (2022). Elucidation of partial activation of cannabinoid receptor type 2 and identification of potential partial agonists: molecular dynamics simulation and structure-based virtual screening. Comput. Biol. Chem.99, 107723. 10.1016/J.COMPBIOLCHEM.2022.107723

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  Valdés-TresancoM. S.Valdés-TresancoM. E.ValienteP. A.MorenoE. (2021). Gmx_MMPBSA: a new tool to perform end-state free energy calculations with GROMACS. J. Chem. Theory Comput.17 (10), 6281–6291. 10.1021/acs.jctc.1c00645

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  VeberD. F.JohnsonS. R.ChengH.-Y.SmithB. R.WardK. W.KoppleK. D. (2002). Molecular properties that influence the oral bioavailability of drug candidates. J. Medicinal Chemistry45 (12), 2615–2623. 10.1021/jm020017n

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  VickeryO. N.StansfeldP. J. (2021). CG2AT2: an enhanced fragment-based approach for serial multi-scale molecular dynamics simulations. J. Chem. Theory Comput.17 (10), 6472–6482. 10.1021/ACS.JCTC.1C00295/ASSET/IMAGES/LARGE/CT1C00295_0008.JPEG

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 60

  WangY.LiuM.GaoJ. (2020). Enhanced receptor binding of SARS-CoV-2 through networks of hydrogen-bonding and hydrophobic interactions. Proc. Natl. Acad. Sci.117 (25), 13967–13974. 10.1073/pnas.2008209117

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  XingC.ZhuangY.XuT. H.FengZ.ZhouX. E.ChenM.et al (2020). Cryo-EM structure of the human cannabinoid receptor CB2-Gi signaling complex. Cell180 (4), 645–654.e613. 10.1016/j.cell.2020.01.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  ZdrazilB.FelixE.HunterF.MannersE. J.BlackshawJ.CorbettS.et al (2024). The ChEMBL database in 2023: a drug discovery platform spanning multiple bioactivity data types and time periods. Nucleic Acids Research52 (D1), D1180–D1192. 10.1093/nar/gkad1004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  ZengC.DubreuilM.LarochelleM. R.LuN.WeiJ.ChoiH. K.et al (2019). Association of tramadol with all-cause mortality among patients with osteoarthritis. JAMA321 (10), 969–982. 10.1001/JAMA.2019.1347

  • Pubmed Abstract
  • CrossRef
  • Google Scholar