Systems Pharmacology Modeling Identifies a Novel Treatment Strategy for Bortezomib-Induced Neuropathic Pain

Chemotherapy-induced peripheral neurotoxicity is a common dose-limiting side effect of several cancer chemotherapeutic agents, and no effective therapies exist. Here we constructed a systems pharmacology model of intracellular signaling in peripheral neurons to identify novel drug targets for preventing peripheral neuropathy associated with proteasome inhibitors. Model predictions suggested the combinatorial inhibition of TNFα, NMDA receptors, and reactive oxygen species should prevent proteasome inhibitor-induced neuronal apoptosis. Dexanabinol, an inhibitor of all three targets, partially restored bortezomib-induced reduction of proximal action potential amplitude and distal nerve conduction velocity in vitro and prevented bortezomib-induced mechanical allodynia and thermal hyperalgesia in rats, including a partial recovery of intraepidermal nerve fiber density. Dexanabinol failed to restore bortezomib-induced decreases in electrophysiological endpoints in rats, and it did not compromise bortezomib anti-cancer effects in U266 multiple myeloma cells and a murine xenograft model. Owing to its favorable safety profile in humans and preclinical efficacy, dexanabinol might represent a treatment option for bortezomib-induced neuropathic pain.


INTRODUCTION
Chemotherapy-induced peripheral neuropathy (CIPN) is a common adverse side effect of cancer chemotherapy that occurs in approximately 30-40% of patients (Wolf et al., 2008;Staff et al., 2017). CIPN typically manifests as a tingling and numbness sensation in the extremities of the body, which can develop into a burning pain (Wolf et al., 2008). Muscle weakness, fatigue, and irreversible nerve damage can be life debilitating for many patients. Improvement in cancer therapy has resulted in a decrease in cancer death rates and an increase in the number of adult and pediatric cancer survivors who may struggle to manage long-term adverse side effects from treatment such as CIPN (Miaskowski et al., 2017;Kandula et al., 2018). An expert panel of the American Society of Clinical Oncology reported that there were no effective agents for preventing CIPN in cancer patients treated with neurotoxic chemotherapy based on a review of clinical trials conducted between 1946 to April 2013 (Hershman et al., 2014). Several nutraceuticals have been evaluated for their potential utility in CIPN, such as vitamin E, goshajinkigan, acetyl-L-carnitine, and alpha-lipoic acid; however, none were effective. Only one drug, duloxetine, has been moderately recommended to treat painful peripheral neuropathy (Hershman et al., 2014). Opioids, cannabis, and gabapentinoids have been used to treat neuropathic pain associated with CIPN. However, these symptomatic therapies exhibit limited efficacy, often associated with an addictive potential and adverse effects.
The molecular mechanisms of neurotoxicity by chemotherapeutics are multifactorial and complex, and there are a wide range of cellular components and processes that are disrupted. Some mechanisms may be exclusive to specific classes of chemotherapeutics (e.g., microtubule or proteasome inhibitors); however, many molecular pathways of toxicity are shared across neurotoxic agents. Some of the possible mechanisms involved in the pathogenesis of CIPN and neuropathic pain include mitochondrial changes, alterations in ion channels and ionic current, inflammation, oxidative and endoplasmic reticulum (ER) stress, activation of the intrinsic apoptosis pathway, MAPK pathway alterations, and the modulation of NMDA receptors (Jaggi and Singh, 2012). Mitochondrial dysfunction has emerged as a fundamental cause of CIPN, which consists of aberrant Ca 2+ signaling, oxidative stress, and neuronal apoptosis. Mitochondria in the sensory axons of dorsal root ganglion (DRG) neurons become swollen and vacuolated when exposed to paclitaxel, oxaliplatin, and bortezomib. Opening of the mitochondrial permeability transition pore (mPTP) increases the permeability of the inner mitochondrial membrane resulting in a loss of membrane potential, mitochondrial swelling, decreased energy production, and rupture of the outer mitochondrial membrane, which initiates apoptosis through the release of pro-apoptotic proteins. Neuroinflammation is another important mechanism governing the development of CIPN (Lees et al., 2017). Upon peripheral nerve injury, resident macrophages, Schwann cells, satellite glial cells, and recruited immune cells secrete prostaglandins, chemokines, and cytokines. Neuroinflammatory mediators of the nerve microenvironment drive the complex interplay between cells of the neuroimmune system to promote adaptive (survival and growth) and maladaptive (neuropathic pain) responses. Several of the highly neurotoxic chemotherapeutics (paclitaxel, oxaliplatin, vincristine, and bortezomib) increase pro-inflammatory cytokines (TNFα and IL-1β) and downregulate anti-inflammatory cytokines (IL-10) in the DRG and spinal cord. Pro-inflammatory cytokines (TNFα, IL-1β, IL-6) and chemokines (CCL2) might be potential biomarkers and targets for predicting/preventing CIPNrelated pain.
The multifactorial nature by which chemotherapeutics induce peripheral neurotoxicity warrants a systems pharmacology approach to evaluate the role of complex physiological processes involved in the development of CIPN and to identify potential treatments. In this study, a network-based systems pharmacology model of signal transduction and gene regulatory processes in peripheral neurons was constructed to predict the intracellular toxicodynamic effects of proteasome inhibitors and to identify drug targets for CIPN. Since the modulation of a single target will likely be inadequate for preventing CIPN, analyses that consider the entire topology and dynamics of the network were performed to identify combinatorial treatment strategies. A novel targeting approach was identified that consists of a TNFα inhibitor, NMDA receptor antagonist, and reactive oxygen species (ROS) inhibitor. Dexanabinol, a synthetic cannabinoid derivative, has been reported to inhibit all three targets (Eshhar et al., 1993;Nadler et al., 1993;Eshhar et al., 1995;Shohami et al., 1997). In this study, the neuroprotective effects of dexanabinol were evaluated in preclinical models of bortezomib-induced neurotoxicity.

Systems Pharmacology Model Development and Qualification
A bottom-up model development approach was utilized (Wang et al., 2012;Bloomingdale et al., 2018), and network construction was initiated using genes specific to drugs known to induce peripheral neuropathy, which were identified from a prior interaction network (Hur et al., 2014). A list of 230 genes was uploaded to DAVID and mapped to pathways in KEGG to identify significant gene-associated pathways. Signaling pathways were cross-referenced with literature to select pathways associated with peripheral neuropathy and not therapeutic indications. Three major pathways were selected (neurotrophin, mitogen-activated protein kinase, and apoptosis) as a foundational base network. Additional overlapping pathways were also included: TNFα, calcium signaling, PI3K-Akt, NFκB, p53, RAS, and protein processing in the endoplasmic reticulum. Other databases (e.g., Reactome and WikiPathways) were also used in model construction. Although peripheral neurons are post-mitotic, pathways related to the cell cycle and proliferation were included since aberrant cell cycle reentry can cause neuronal death.
Model qualification was achieved by comparing model predictions with an external transcriptomic dataset not used during network construction: Gene Expression Omnibus (GEO), GSE10470, in which gene expression was measured in neurons isolated from ATF4 knockout mice with and without exposure to homocysteate, an NMDA receptor agonist. Differentially expressed genes for each experimental condition were identified using Bioconductor packages (Biobase v2.30.0,GEOquery v2.40.0,limma v3.26.8). Data were log transformed using a base of two. Linear least squares regression was performed for each gene, and differentially expressed genes were identified using an empirical Bayes method with Benjamini and Hochberg FDR adjusted p-values. Network predictions were performed using normalized HillCube functions in Odefy, a MATLABbased toolbox to convert Boolean models to a series of ordinary differential equations (Krumsiek et al., 2010). ATF4 and NMDA were fixed to zero and one to represent the gene knockout and presence of homocysteate. Simulations were performed until steady-state was reached. Network predictions for the change in node expression were reported as increasing when the initial state changes from 0 or 1 to a final state of 1, and decreasing when the initial state changes from 0 or 1 to a final state of 0. Network predictions were compared to the Log2 foldchange in differentially expressed genes. In cases where multiple probes exist for the same gene, the probe with the smallest adjusted p-value was selected for comparison. Averaging gene expression across multiple probes provided identical results.

Network Model Analyses
Bortezomib was assumed to only inhibit the proteasome. Boolean logic functions were converted to normalized HillCube differential equations using Odefy, and default parameter values were used. Network simulations were performed for 25time steps (arbitrary units), which was sufficient to reach steadystate.
Minimal intervention analysis (MIA) was performed using CellNetAnalyzer, a MATLAB-based toolbox. The intervention goal was to identify two-and three-node interventions that retain apoptosis in an OFF (zero) state in the presence of proteasome inhibition (proteasome fixed to zero). Nerve growth factor and brain derived neurotrophin factors were also assumed to be present (fixed to one).
An attractor analysis using a synchronous updating scheme and 10 6 initial start states was performed for eight network perturbations. Proteasome was fixed to zero for all eight analyses to represent the effect of bortezomib. Relative activation frequencies were determined for all network components and apoptosis, which was calculated as the sum of node activation states in the attractors of each simulation divided by the total number of simulations performed. All analyses were performed using BoolNet.

Neuronal and Multiple Myeloma Cell Culture
Neuronal SH-SY5Y cells were grown at 37°C and 5% CO 2 in DMEM:F12 media supplemented with 10% fetal bovine serum (10,000 U/mL) and 1% penicillin-streptomycin (10,000 μg/ml) until a sufficient number of cells were obtained. SH-SY5Y cells were seeded at a density of 100,000 cells/well in 96 well plates and allowed to equilibrate for 24 h. Bortezomib (0.001, 0.01, 0.1, 1, 10, 100, 1,000 nM) and dexanabinol (1 and 10 μM) were added to their respective wells. Concentration ranges were selected in order to achieve a wide range around the IC 50 of bortezomib and a potential therapeutic range of dexanabinol. Cell viability was measured using a colorimetric reagent (WST-1) (Berridge et al., 2005). After 48 h of drug exposure, 10% v/v of WST-1 was added to each well. Plates were briefly shaken and placed into an incubator for 2 h. Absorbance was measured using a SpectraMax 190 microplate reader at a wavelength of 450 nm (690 nm reference wavelength).
U266 myeloma cells were grown at 37°C and 5% CO 2 in RPMI-1640 media supplemented with 15% fetal bovine serum (10,000 U/mL) and 1% penicillin-streptomycin (10,000 μg/mL) until a sufficient number of cells were obtained. U266 cells were seeded at a density of 25,000 cells/well in 96 well plates and were allowed to equilibrate for 24 h. Bortezomib (0.001, 0.01, 0.1, 1, 10, 100, 1,000 nM) and dexanabinol (1, 3, 10, 30, 100, 300 μM) were added to their respective wells, and concentration ranges were selected to achieve a wide range around the IC 50 of each drug. Cell viability was measured after 24, 48, or 72 h of drug exposure as described above.

Pharmacodynamic Modeling of Bortezomib and Dexanabinol In Vitro Cytotoxicity
For single agents, the decrease in cell viability (%) as a function of drug concentration was modeled using an inhibitory Hill function (Black and Leff, 1983). For the combinatorial effect of two drugs, cell viability was modeled using a modified version of the Ariens equation (Chakraborty and Jusko, 2002). An interaction parameter (Ψ vit ) was incorporated on the bortezomib IC 50 to account for synergistic (Ψ vit < 1), antagonistic (Ψ vit > 1), or additive (Ψ vit 1) effects on U266 cell viability. Parameters were estimated from single agent treatments and fixed while modeling interaction data to estimate Ψ vit .

Dexanabinol Effects on LPS-Induced TNFα Production in Macrophages
WEHI-13VAR fibroblast cells, a TNFα sensitive cell line derived from a mouse fibrosarcoma, were grown in RPMI-1640 culture medium containing 2 mM L-glutamine, 10% fetal bovine serum, and 3 μg/ml gentamicin in T75 flasks at 37°C, 95% relative humidity, and 5% CO 2 . Cells were cultured to approximately 90% confluency and were always below passage 25 to avoid loss of TNFα sensitivity. Cells were prepared by detaching with 0.25% trypsin and 0.02% EDTA and resuspending in culture medium supplemented with 1 μg/ml actinomycin D to a concentration of 500,000 cells/ml. 100 µl of cell suspension were added to each well of a flat-bottom 96well plate containing 100 µl of 2-fold serial dilutions of unknown samples, in duplicate, or known concentrations of rat recombinant TNFα standards in diluting medium (RPMI-1640, 2 mM L-glutamine, 1% fetal bovine serum, and 15 mM HEPES). Following 20 h of incubation, 10 µl of Cell Proliferation Reagent WST-1 (Roche Diagnostics, Indianapolis, IN, United States) in diluting medium was added to each well. WST-1 counting solution was used as a cell viability indicator, which was quantified spectrophotometrically (Berridge et al., 2005). After 4 h of incubation, absorbance at 440/700 nm was measured using a SpectraMax 96 microplate reader with SoftMax Pro v.4.0 acquisition and analysis software (MDS Analytical Technologies, Sunnyvale, CA, United States). A standard curve for absorbance (OD440-OD700) vs TNFα was plotted. The assay detection limit is approximately 1 pg/ml (Eskandari et al., 1990).

Microphysiological Peripheral Nerve Model
All animal handling and tissue harvesting procedures were performed according to guidelines set by NIH and the Institutional Animal Care and Use Committee (IACUC) at Tulane University. DRG from an E15 rat litter were harvested, pooled, dissociated in trypsin, and passed through a 40 µM nylon mesh filter. Cells were then seeded in a 96 well ultra-low attachment treated spheroid microplate (Corning Inc., Corning, NY, United States) at a concentration of 45,000 cells per well in spheroid formation media composed of Neurobasal Medium supplemented with 2% v/v B27 supplement, 1% v/v GlutaMAX, 20 ng/ml nerve growth factor 2.5 S native mouse protein (NGF), and 1% v/v antibiotic/antimycotic solution (all from Thermo-Fisher Scientific, Waltham, MA, United States). Microplates were centrifuged at 500 g for 5 min and spheroids aggregate after overnight incubation at 37°C and 5% CO 2 .
Nerve impulse conduction was evaluated using extracellular field potential recording. All recordings were performed in continuously perfused, room temperature, bicarbonatebuffered, artificial cerebrospinal fluid (ACSF) bubbled with 95% O 2 , 5% CO 2 as described previously (Huval et al., 2015). Electrical stimulation of nerve cultures and recording of compound action potential (CAP) propagation was performed similar to a previous description (Sharma et al., 2019). CAP propagation into the neuronal cell body cluster was recorded after nerite stimulation at two locations, 1.5 and 3 mm distal from the recording site. For each location, a stimulation pulse height of 15 V and pulse width of 200 μs was replicated ten times at 0.5 Hz and the ten resulting traces were averaged for analysis. Action potential amplitude (APA) is measured as the most negative potential recorded in the CAP trace. Recording distance was defined as the distance between the stimulating and recording electrodes. Latency to CAP onset was measured as the time elapsed between stimulation and the beginning of the CAP signal in the cell body cluster. Nerve conduction velocity (NCV) was defined as the recording distance divided by the latency to CAP onset.

Bortezomib-Induced Peripheral Neuropathy Rat Model
A well-established bortezomib-induced peripheral neuropathy rat model was used to evaluate the neuroprotective effects of dexanabinol (Meregalli et al., 2010). Thirty younger adult female Wistar rats, (175-200 g; Envigo, Udine, Italy) randomized in 10 animals/group were used. Care and husbandry of animals were in conformity with the institutional guidelines in compliance with national (D.L. n. 26/2014) and international laws (EEC Council Directive 86/609, OJ L358, 1, Dec. 12, 1987; Guide for the Care and Use of Laboratory Animals, U.S. National Research Council, 1996). The experiments were examined and approved by the Ethics Committee of the University of Milano-Bicocca and Ministry of Health (approval number: 0007583/19 and 842/ 2018-PR, respectively). All rats were housed in a limited-access animal facility at the University of Milano-Bicocca, where room temperature, humidity, and artificial lighting were set to previously determined conditions. A veterinarian with specific expertise in animal studies and the authority to withdraw the animal from the experiment examined any animal showing signs of distress. Bortezomib and dexanabinol were purchased from LC Laboratories (Woburn, MA, United States, Canada) and Cayman Chemical Company (Ann Arbor, MI, United States).
Animals were allocated to experimental groups based on their baseline responses to the behavioral and neurophysiological studies, randomizing them into three groups containing 10 rats each: 1) dexanabinol vehicle-treated controls (vehicle) 2) bortezomib (Bort), and 3) dexanabinol + bortezomib (Dex + Bort). Dexanabinol was administered 30 min before bortezomibinjections at 9:00 am, and bortezomib or vehicle-treatments were performed between 9:30 to 11 am.
Dynamic and plantar aesthesiometry were performed under blinded conditions and were carried out by different researchers (C.M and L.M, respectively) between 9:00 to 12:00 am, as well as Frontiers in Pharmacology | www.frontiersin.org January 2022 | Volume 12 | Article 817236 intraepidermal unmyelinated axons were counted in a blinded fashion (A.Ch). Fresh solutions were prepared before each administration. Bortezomib was dissolved in a solution of 10% Tween 80, 10% ethanol 100%, and 80% saline solution and was administered intravenously (IV) (0.2 mg/kg) three times a week for 8 weeks. Dexanabinol was formulated as a cremophor/ethanol solution (70%/30%) after a 1:20 dilution in saline. Dexanabinol was administered intraperitoneally (IP) (10 mg/kg) approximately 30 min prior to the administration of bortezomib. The general condition of the animals was assessed daily during the treatment period, and body weights were recorded twice weekly for toxicity assessment.
To assess the development of neurotoxicity sensory NCV and SAP amplitude were measured in caudal nerve using an electromyography apparatus (Myto2 ABN Neuro, Firenze, Italy) according to previous study (Monza et al., 2021). Briefly, distal caudal nerve was recorded orthodromically: recording cathode and anode were placed at 6 and 5 cm from the tip of the tail, the ground electrode at 2.5 cm from it, and stimulating anode and cathode respectively at 2 and 1 cm. Filters were set between 20 Hz and 3 KHz for sensory recordings. Animals were kept under deep isoflurane anesthesia, and body temperatures were kept constant at 34.5 ± 0.5°C for the whole procedure.
Pain sensitivity as a behavioral outcome was assessed prior to treatment (baseline) and after four and 8 weeks of BTZadministration using the dynamic plantar aesthensiometer test (model 37450; Ugo Basile Biological Instruments, Comerio, Italy) and the dynamic plantar analgesiometer (model 37370; Ugo Basile Biological Instruments) (Chiorazzi et al., 2018). Mechanical allodynia was evaluated by quantifying the mechanical-escape behavior for 20 s after a metal filament was applied to the plantar skin of the hind paw according to methods as previously described (Chiorazzi et al., 2018). Mechanical threshold was automatically assessed by the mean value of six repeated application, measuring the paw withdrawal threshold (as withdrawal threshold expressed in grams) against stimulation with von Frey filaments, which exerted a progressively increasing puncture pressure reaching up to 50 g. Two hours after dynamic test evaluation, the response to noxious thermal stimulus was determined by a movable infrared radiant heat source (intensity of 40 IR), which was placed directly under the plantar surface on the hind paw and time taken for hind paw withdrawal was monitored (as withdrawal latency expressed in seconds). The nociceptive threshold response to thermal stimulus was determined by calculating the mean value of four repeated measurements. For both measures, a cut-off time of 25 s was used to prevent tissue damage in the absence of a response. All tests were operated by an experimenter blinded to groupings and in the controlled behavioral test room.
The damage of peripheral small nerve fibers, an analysis of intraepidermal nerve fiber (IENF) density, was performed at the end of the 8 weeks of drug treatment as previously described (Canta et al., 2016). Briefly, intraepidermal unmyelinated axons were counted in a blinded fashion, and the IENF density was obtained by immunochemistry analysis counting polyclonal antiprotein gene product 9.5 (PGP 9.5; GeneTex, Irvine, CA, United States)-positive fibers IENF/length of the epidermis (mm) (Canta et al., 2016).

Multiple Myeloma Mouse Model
All animal protocols were approved by the Institutional Animal Care and Use Committee (IACUC) at Roswell Park Comprehensive Cancer Center. SCID mice were obtained from Laboratory Animal Resource of Roswell Park Comprehensive Cancer Center. Human multiple myeloma MM1.S cell line (gift from Dr Stephen Rosen, Chicago, IL, United States) was maintained in RPMI 1640 media containing 10% fetal bovine serum at 37°C under 5% CO 2 . Twenty SCID female mice at age of 4 weeks were irradiated with 300 Rad using a Mark II Cesium irradiator. Twenty-4 hours after irradiation, five million MM1.S cells were injected subcutaneously in the left ventral flank. Animals were monitored for tumor formation, and tumor sizes were measured with a digital caliper twice weekly. Tumor volumes were calculated using the modified ellipsoidal formula: tumor volume 0.5 (length × width × width). When tumor volumes reached approximately 100 mm 3 , animals were divided to four groups with similar average tumor volumes. Treatments were administered twice a week via IP injection at the doses of 1 and 10 mg/kg for bortezomib and dexanabinol. The study was conducted for 4 weeks and mice were sacrificed when tumors grew to the limit enforced by IACUC.

Pharmacokinetic and Pharmacodynamic Modeling
A PBPK model of bortezomib in mice (Zhang and Mager, 2015) was used to simulate a multiple dosing regimen of 1 mg/kg bortezomib administered IV to mice twice weekly. The average body weight of all mice in the bortezomib treatment group was used to determine the dose for each administration. Dexanabinol PK in mice has not been reported; therefore, dexanabinol PK in rats and humans (Brewster et al., 1997;Gershkovich et al., 2007) was allometrically scaled to predict mouse dexanabinol PK. A two-compartment PK model was fit to rat and human plasma concentrations of dexanabinol, and a power law relationship was used for parameter scaling: with α and β calculated for each PK model parameter (X) using species body weight (BW) (Boxenbaum, 1982). The values for α and β, along with the average body weight of mice in the dexanabinol treated group, were used to predict plasma dexanabinol PK in mice.
Bortezomib and dexanabinol PK were fixed when estimating pharmacodynamic model parameters. The rate of change of tumor volume (TV) was defined by the following differential equations: Frontiers in Pharmacology | www.frontiersin.org January 2022 | Volume 12 | Article 817236 Tumor volume in the absence of drug was modeled using an exponential growth function (Norton and Simon, 1977). Drug effect was modeled as second-order elimination processes governed by the rate constants K Bort and K Dex . A function that decreases the rate of elimination of tumor volume by bortezomib was included (k res ), which improved model performance. The interaction term (Ψ viv ) was incorporated on the parameter governing bortezomib resistance (k res ). Ψ viv was introduced to account for synergistic (Ψ viv < 1), antagonistic (Ψ viv > 1), or additive (Ψ viv 1) pharmacodynamic interactions. PK/PD modeling was performed using ADAPT 5 (D' Argenio et al., 2009).

Statistical Analysis
Statistical analyses were performed using GraphPad Prism version 7.04. The change in nerve-on-a-chip NCV and APA for various concentrations of bortezomib was assessed using oneway ANOVA with Dunnet's correction for multiple comparisons. Comparisons were made for each group with respect to the vehicle control. The change in nerve-on-a-chip NCV and APA for the drug combination was assessed using one-way ANOVA with Tukey's correction for multiple comparisons. Drug pharmacodynamics on tumor volume in the multiple myeloma SCID mouse model relative to control was tested via multiple t-tests with Holm-Sidak Bonferroni correction (comparisons were made across all groups). IENF density, pain thresholds, and neurophysiological assessments in the rat model were analyzed using a Kruskal-Wallis test with Dunn's correction for multiple comparisons across all groups.

Systems Pharmacology Model Development and Qualification
A Boolean model of neuronal signaling was developed, which contains 131 nodes/components and 252 edges or reactions (Supplementary Figure S1). All model equations, initial conditions, and supporting references are provided in Supplementary Table S1. The model contains nine input nodes (BDNF, NGF, IL-1β, TFGβ, TNFα, IL-6, FasL, NMDA, and proteasome), which primarily consist of ligands that bind to 9 cell membrane receptors (NGFR, TrkA, TrkB, IL-1R, TFGβR, TNFαR, IL-6R, FasR, and NMDAR). Model inputs are connected to four endpoints (nerve damage, microtubule transport, apoptosis, and proliferation) through a complex intracellular network.
Network qualification was performed using a gene expression profile from an external microarray dataset of neurons isolated from ATF4 knockout mice with and without exposure to homocysteate, an NMDA receptor agonist. ATF4 is a central node in our network model and appears to be important component in oxidative/ER stress pathways. The ATF4 knockout provides a reasonable indirect dataset to test the connectivity of the final network model. ATF4 is downstream of NMDAR and is involved in the ROS signaling pathway, and the ATF4 dataset was considered suitable for model qualification purposes. Network-predicted changes in node activation states were compared to the Log 2 fold-change in differentially expressed genes. There was a total of 2,113 differentially expressed genes between ATF4 knockout murine neurons exposed to an NMDA agonist and treatment naïve wild-type neurons. Nineteen genes that exhibited differential expression overlapped with nodes in the network model. Network predictions for the change in node expression were recorded as "increasing" when the initial state changed from 0 or 1 to a final state of 1 and "decreasing" when the initial state changed from 0 or 1 to a final state of 0. Qualitative model predictions agreed with the directionality of experimental observations for 16 out of 19 nodes (84%) (Supplementary Figure S2).

Systems Model-Predicted Treatment Strategies for CIPN
Two types of network analyses, minimal intervention analysis (MIA) and attractor analysis, were performed to identify treatment strategies for proteasome-inhibitor induced peripheral neuropathy. MIA identified several targets/ pathways, whereas attractor analysis enabled the prioritization of intervention targets based upon their relative importance.
MIA identified 224 potential combinatorial targets that might prevent neuronal apoptosis induced by proteasome inhibition. The potential strategies consisted of 109 two-target and 115 three-target interventions. The frequency of node interventions relative to the total number of intervention sets (up to two target combinations) are shown in Supplementary Figure S3A, and all two-target intervention sets are shown in Supplementary Figure S3B. Components of the TNFα signaling pathway, (i.e., TNFα, TNFαR, TRADD, and FADD) were the most frequently identified species. A TNFα inhibitor in combination with a pharmacological agent that modulates one of the 27 targets listed under Target 2 in Supplementary Figure  S3B is predicted to prevent proteasome inhibitor induced neuronal apoptosis. A TNFα inhibitor in combination with an inhibitor of ROS was identified as a two-target intervention set. The frequency of node interventions relative to the total number of intervention sets (up to three combinations) are shown in Supplementary Figure S3C. All three-target intervention sets are shown in Supplementary Figure S3D. Components of the TNFα signaling pathway were also prominent across the three-target intervention sets. At least one node in the TNFα signaling pathway (TNFα, TNFαR, TRADD, and FADD) was identified in 88 of the 115 threetarget intervention sets. Three-target intervention sets are classified into five distinct groups (Supplementary Figure  S3D), and the first group in particular consists of the inhibition of nodes in the TNFα signaling pathway, NMDA pathway (NMDA, NMDAR, NO, nNos), and MAPK pathway (MEKK1, MKK4, MKK7, TAK1).
Network simulations were performed to assess the intraneuronal toxicodynamics of proteasome inhibition. Predicted normalized expression profiles of select cellular components is displayed in Figure 1A. Model predictions for all cellular components included in the network are shown in Supplementary Figure S4. Proteasome inhibition resulted in the production of ROS, which activated ER stress and the UPR. Several transcription factors increased in expression (e.g., ATF4 and CHOP), pro-survival proteins (Bcl-2, Bcl-xL, XIAP) decreased, and pro-apoptotic proteins (Bad, Bax, Cyc, Casp3, and Casp9) increased in expression. Lastly, proteasome inhibition resulted in the activation of neuronal apoptosis.
Based on TNFα and NMDA pathways being prominently featured in the MIA (Supplementary Figure S3), and simulations suggesting a role for ROS ( Figure 1A), an attractor analysis was performed to predict the role of a TNFα inhibitor, NMDA receptor antagonist, and ROS inhibitor, alone and in combination, in decreasing the frequency of neuronal apoptosis ( Figure 1B). For each treatment group, 10 6 model simulations were performed, and the relative activation frequencies of network nodes are depicted as a heatmap. In the first column (bortezomib alone), nodes in the TNFα and NMDA pathways were equal to 0.5 and ROS was equal to 1.
Hence, the TNFα/NMDA pathways are active at steady-state in half of the simulations, and the oxidative stress pathway is active at steady-state throughout all simulations. In columns 2-4, the relative activation frequencies of nodes in the NMDA, ROS, and TNFα pathways decrease to zero in the presence of their respective inhibitors. A trend of decreasing pro-apoptotic (Bax, Casp7, Casp3, Cyc, Casp9) and increasing pro-survival (Bcl-2, XIAP, Bcl-xL) species was observed. This trend becomes more pronounced for two-target combinations (columns 5-7) and is highly evident in the three-target combination (column 8). The relative activation frequency of neuronal apoptosis (ν apoptosis ) was determined for the eight treatment arms ( Figure 1C). ν apoptosis was 0.69 in the bortezomib only group, which decreased to 0.63, 0.51, and 0.38 in the presence of single-agent treatments with an NMDA inhibitor, ROS inhibitor, or TNFα inhibitor. ν apoptosis decreased further in the two-target combinations and was completely inhibited in the presence of all three inhibitors ( Figure 1C).

Assessment of Bortezomib and Dexanabinol In Vitro Cytotoxicity
SH-SY5Y neuroblastoma cells were used to assess the neuroprotective effect of dexanabinol on preventing bortezomib-induced cytotoxicity. The viability of SH-SY5Y FIGURE 1 | Network simulations to identify combinatorial treatment strategies for BIPN. (A) Network model simulations of select intracellular components in the presence of proteasome inhibition. The Boolean network model was converted to normalized HillCube differential equations using Odefy, and default parameter values were used (tau 1; k 0.5; n 3). Simulations were performed for 25-time steps (arbitrary units). Results from attractor analysis for the relative activation frequency of (B) network nodes and (C) neuronal apoptosis across identified attractors for eight network perturbations: proteasome inhibition (bortezomib) in combination with a TNFα inhibitor, NMDA receptor antagonist, and/or ROS inhibitor. Monotherapies, two-target combinations, and the three-target combinations are highlighted in red, blue, and purple.
Frontiers in Pharmacology | www.frontiersin.org January 2022 | Volume 12 | Article 817236 cells at 48 h was determined for a range of bortezomib concentrations in the absence and presence of dexanabinol (1 and 10 μM; Figure 2A). A neuroprotective effect was not observed in SH-SY5Y cells, with the IC 50 of bortezomib alone (4.40 nM) being comparable to the IC 50 estimates in the presence of 1 and 10 μM of dexanabinol (7.00 and 4.62 nM). Greater concentrations of dexanabinol (≥25 μM) were shown to be cytotoxic (data not shown). U266 multiple myeloma cells were used to assess the potential anti-cancer effects of dexanabinol and whether it would prevent bortezomib-induced cytotoxicity. The viability of U266 multiple myeloma cells was determined for a range of bortezomib and dexanabinol concentrations. A BLISS independence model was used to assess the nature of the pharmacological interaction between the two compounds at 24 ( Figure 2B Tables S2-S4). An antagonistic relationship was observed around the IC 50 of bortezomib for low concentrations of dexanabinol (≤10 μM) ( Figure 2B). This antagonistic relationship appeared to be more pronounced at 24 and 48 h compared to 72 h. The timedependent shift from an antagonistic to additive pharmacological relationship was quantified through a decrease in Ψ vit from 2.71 (24 h) to 2.10 (48 h) to 1.40 (72 h) ( Supplementary Tables S2-S4).

Dexanabinol Inhibits Macrophage TNFα Production in a Biphasic Manner
WEHI-13VAR fibroblast cells were used to analyze supernatants from primary rat macrophage cultures for the presence of biologically active TNFα (Khabar et al., 1995;Ignatowski et al., 2000). Macrophages were stimulated with LPS to produce TNFα, which was inhibited in the presence of dexanabinol ( Figure 2C). The concentration-effect relationship was bell-shaped, and maximal inhibition of TNFα occurred at 1 μM.

Dexanabinol Exhibits Neuroprotective Effects at Low Concentrations Using Nerve-On-a-Chip
The efficacy of dexanabinol in preventing the neurotoxic effects of bortezomib was evaluated using a microphysiological model of a rat peripheral nerve, referred to as nerve-on-a-chip (Huval et al., 2015). Nerve conduction velocity (NCV) and action potential amplitude (APA) was measured for various concentrations of bortezomib to obtain a concentration-effect relationship (Supplementary Figure S5). Proximal and distal recordings of NCV and APA exhibited a decreasing trend for increasing concentrations of bortezomib. Proximal NCV at 1,000 nM of bortezomib decreased; however, this was not statistically significant (p 0.054; Supplementary Figure S5C). There was FIGURE 2 | Assessment of in vitro cytotoxicity of bortezomib and dexanabinol. Cell viability was measured using WST-1. (A) Cytotoxicity of SH-SY5Y cells exposed to various concentrations of bortezomib in the absence (red) or presence of 1 μM (green) and 10 μM (blue) dexanabinol. (B) Cytotoxicity of U266 multiple myeloma cells exposed to a range of bortezomib and dexanabinol concentrations for (left) 24, (middle) 48, and (right) 72 h. Antagonistic, additive, and synergistic pharmacological relationships are colored in red, green, and blue. (C) Concentration-effect relationship of dexanabinol on LPS-induced TNFα production in primary rat macrophages.
Frontiers in Pharmacology | www.frontiersin.org January 2022 | Volume 12 | Article 817236 a significant decrease in distal NCV after 48 h of bortezomib exposure at 100 and 1,000 nM (p < 0.05; Supplementary Figure S5D). The effects of dexanabinol on bortezomib-induced changes in APA ( Figures 3A,B) and NCV ( Figures 3C,D) were evaluated. Neurophysiological endpoints were normalized to their control values to account for batch-to-batch variability. The effects of 10 and 25 μM dexanabinol on NCV and APA were not significantly different from control, which indicates that dexanabinol is not neurotoxic at the selected concentrations. Nerves exposed to 100 nM of bortezomib for 48 h exhibited a significant decrease in the proximal amplitude and distal NCV (Figures 3A,D). The treatment group of 100 nM of bortezomib in combination with 10 μM of dexanabinol resulted in NCVs and amplitudes that were no longer significantly different than control. However, there was no significant difference observed in distal NCV between bortezomib alone and in combination with 10 μM of dexanabinol (p value 0.070). The combination of bortezomib with a greater concentration of dexanabinol (25 μM) trended towards a further decrease in APA and NCV, which indicates a possible enhanced neurotoxic effect at greater dexanabinol concentrations. The biphasic inhibition of TNFα by dexanabinol ( Figure 2C) might partially explain the improvement in proximal APA and distal NCV at lower dexanabinol exposure (10 μM) but not for the higher concentration (25 μM).

Dexanabinol Prevents Bortezomib-Induced Allodynia and Thermal Hyperalgesia in Rats
The onset of peripheral neuropathy manifested after 4 weeks of chemotherapy with a significant decrease in distal, but not proximal, sensory action potential (SAP) amplitude (p < 0.001) and sensory NCV (p < 0.01) in bortezomib-treated groups compared to vehicle (Figure 4). Dexanabinol treatment did not prevent bortezomib-induced decreases in neurophysiological endpoints (Figure 4). After 8 weeks of bortezomib therapy, there were significant decreases in proximal caudal SAP amplitude (p < 0.01), distal caudal SAP amplitude (p < 0.001), and distal caudal sensory NCV (p < 0.01) ( Figures 4A,B,D). The combination of dexanabinol and bortezomib resulted in further decreases in proximal caudal SAP amplitude than bortezomib alone with respect to vehicle ( Figure 4A) (p < 0.001 vs p < 0.01). Similarly, co-treated groups with Dex + Bort showed a reduction of both distal caudal SAP (p < 0.001) and sensory NCV (p < 0.01) compared to vehicle ( Figures 4B,D). After four bortezomib-treatments, rats developed mechanical allodynia and thermal hyperalgesia, with a significant decrease in mechanical (p < 0.001; Figure 5A) and thermal (p < 0.01; Figure 5B) thresholds compared to vehicle. In addition, significant differences between bortezomib and the combination of bortezomib and dexanabinol were observed for mechanical (p < 0.05; Figure 5A) and thermal (p < 0.05; Figure 5B) thresholds. After 8 weeks of bortezomib therapy, significant decreases in mechanical (p < 0.01; Figure 5A) and thermal (p < 0.01; Figure 5B) thresholds compared to vehicle were observed. Dexanabinol therapy prevented bortezomib-induced mechanical allodynia and thermal hyperalgesia ( Figures  5A,B). At 8 weeks, significant differences were observed between bortezomib and the combination of bortezomib and dexanabinol for mechanical (p < 0.01; Figure 5A) and thermal (p < 0.01; Figure 5B) thresholds.
Bortezomib exposure for 8 weeks resulted in a significant decrease in intraepidermal nerve fibers (IENF) density (p < 0.001; Figure 6). Treatment with dexanabinol showed a trend toward improving the bortezomib-induced decrease of IENF density, although these changes were not statistically significant ( Figure 6).

Dexanabinol Does Not Interfere With the Anti-Cancer Activity of Bortezomib in a Multiple Myeloma Mouse Model
The pharmacodynamics of dexanabinol on the anti-cancer activity of bortezomib was investigated in a SCID MM1S multiple myeloma mouse model. Pharmacokinetic models were developed to predict the exposure of bortezomib and dexanabinol in mice, which were used to drive the pharmacodynamic effects on tumor volume. A previously developed physiologically-based pharmacokinetic (PBPK) model was used to predict mouse plasma concentrations of bortezomib (1 mg/kg) administered intravenously (IV) twice weekly ( Figure 7A; Supplementary Table S5) (Zhang and Mager, 2015). The predicted PK profile exhibits a rapid distributive phase followed by a prolonged elimination phase. The predicted C max after the last dose was 2.06 μM, and the trough concentration prior to the last dose was 0.083 μM. Dexanabinol PK in mice was scaled down using a model established using literature data from humans and rats. A twocompartment PK model was fit to plasma concentrations of dexanabinol in rats and humans (Supplementary Figure S6). Estimated PK parameters and the respective body weights of rats and humans were used to calculate allometric parameters (α and β), which were used to calculate PK parameters in mice (Supplementary Table S6). The predicted mouse plasma concentration-time profile of dexanabinol (10 mg/kg) Frontiers in Pharmacology | www.frontiersin.org January 2022 | Volume 12 | Article 817236 administered intraperitoneally (IP) twice weekly is shown in Figure 7B. The predicted C max after the last dose was 28.1 μM, and the model predicted half-lives of dexanabinol in humans, rats, and mice of 4.4, 4.1, and 10.8 h (Supplementary Table S6). The pharmacodynamics of dexanabinol and bortezomib on tumor volume is shown in Figure 7C, and model parameter estimates are listed in Supplementary Table S7. Tumor volumes in the dexanabinol alone group were similar to control. Exposures to bortezomib alone and in combination with dexanabinol both resulted in statistically significant decreases in tumor volume. Dexanabinol did not antagonize the anti-cancer effects of bortezomib, and the combination trended towards improved tumor growth inhibition. The performance of the pharmacodynamic model improved after incorporating an exponential decay of the rate constant (K bort ) governing the elimination of tumor volume by bortezomib exposure, suggesting acquired resistance. The exponential decay parameter (k res ) was estimated to be 0.0076 h −1 . An interaction term (Ψ viv ) was placed on k res and was estimated to be 0.883 with a 95% confidence interval of 0.80-0.97, which indicates a synergistic interaction by dexanabinol for slowing the decay rate of tumor volume elimination by bortezomib. Simulations using a prior signaling network model of multiple myeloma (Ramakrishnan and Mager, 2018) also qualitatively supports that dexanabinol does not impair the anti-cancer effect of bortezomib as tumor-associated apoptosis is predicted to remain active with the combination (Supplementary Figure S7).

DISCUSSION
CIPN remains an unmet medical need of increasing importance owing to the lack of available treatment options, the rising number of cancer survivors managing long-term adverse side effects of chemotherapy, and the need of non-opioid analgesics for the treatment of chronic pain (Manchikanti et al., 2012;Miaskowski et al., 2017;Staff et al., 2017). In this study, a systems pharmacology modeling approach enabled the identification of a potential treatment strategy for neuropathic pain associated with proteasome inhibitors, which was tested and confirmed in preclinical models of bortezomib-induced neurotoxicity.
A network-based systems pharmacology model (Supplementary Figure S1) provided a platform for exploring the molecular mechanisms involved in the pathogenesis of CIPN in silico, which identified potential treatment targets. TNFα was identified as an important signaling pathway that must remain inactive to prevent neuronal apoptosis. This supports the emerging concept of neuroinflammation as a major mechanism underlying the development of CIPN (Lees et al., 2017). Network analyses suggested that the combinatorial inhibition of TNFα, NMDAR, and ROS could prevent the activation of proteasome inhibitor-induced neuronal apoptosis. However, FIGURE 5 | Efficacy of dexanabinol for preventing bortezomib-induced mechanical allodynia and thermal hyperalgesia in rats. Mechanical threshold (grams) was determined using dynamic test at baseline, 4 weeks of treatment, and 8 weeks of treatment (A). The withdrawal latency to an infrared heat stimulus (seconds) was determined using a dynamic plantar analgesiometer at baseline, 4 weeks of treatment, and 8 weeks of treatment (B). Bortezomib (Bort) was administered IV (0.2 mg/kg) three times a week, and dexanabinol (Dex) was administered IP (10 mg/kg) three times a week approximately 30 min prior to bortezomib. GraphPad prism v7.04 was used to perform a Kruskal-Wallis test with Dunn's correction for multiple comparisons across all groups. p values are reported as: Bort vs vehicle (Vehicle), ** (<0.01), *** (<0.001) and Bort vs Dex + Bort, ○ (<0.05), ○ ○ (<0.01). there could be potential clinical benefit for the modulation of single or dual targets. Anti-TNFα antibodies have shown to prevent allodynia in rats and partially restore sensory nerve APA in mice treated with bortezomib (Chiorazzi et al., 2013;Ale et al., 2014). In addition, bortezomib has been shown to induce neuropathic pain through the activation of presynaptic NMDA receptors in the spinal cord (Xie et al., 2017). Over the past decade there has been an increased interest in the use of ketamine, an NMDA receptor antagonist, for the treatment of neuropathic pain (Kamp et al., 2019). Several antioxidant therapies have also been investigated clinically for the treatment of CIPN; however, none of these agents have been effective as single agents (Wolf et al., 2008;Hershman et al., 2014). Duloxetine is another polypharmacological agent that has shown partial effectiveness in the treatment of CIPN (Thor and Katofiasc, 1995;Meng et al., 2019). The exact mechanisms of duloxetine neuroprotective effects are unclear, but the drug modulates alpha-1-adrenergic receptors, norepinephrine, serotonin, p38, and NF-kB. Although our final network model does not include all known or potential duloxetine mechanisms of action, the inhibition of p38/NF-kB were identified as potential targets (Supplementary Figure S3), which provides supporting evidence that the neuroprotective properties of duloxetine may be acting through this pathway (Meng et al., 2019). The combination of multiple therapeutics, such as duloxetine, anti-TNFα antibodies, ketamine, and antioxidants could be evaluated to treat CIPN. Alternatively, the generation of optimized polypharmacological drugs utilizing machine learning and mechanism-based modeling could enable the development of novel drugs for CIPN.
The lack of a clinically meaningful therapy for CIPN might result from the focus on single-target interventions rather than accounting for multiple etiological mechanisms. A cocktail of three single-target therapeutics would not likely be clinically viable; therefore, a single drug candidate with a polypharmacological profile was sought that matched the model-predicted combination strategy. Dexanabinol, a cannabinoid derivative, was selected for the treatment of BIPN based on its ability to inhibit TNFα, NMDAR, and ROS (Eshhar et al., 1995;Shohami et al., 1997). The drug also binds and affects several other (off) targets, and the in vitro pharmacodynamic potencies (EC 50 and IC 50 values) of dexanabinol for these targets are summarized in Supplementary Table S8 (Bloomingdale, 2018). Typical plasma dexanabinol concentrations in preclinical species and humans are greater than the reported in vitro potencies (Figure 7 and Supplementary Figure S6). Interest in targeting the endocannabinoid system for the treatment of CIPN has increased owing to preclinical and clinical efficacy in other sensory neuropathies (Abrams et al., 2007;Selvarajah et al., 2010;Blanton et al., 2019). Dexanabinol does not exhibit psychoactive effects and is not prone to cannabinoid receptor desensitization since it does not bind to cannabinoid receptors (Shohami et al., 1997). Dexanabinol was previously investigated for traumatic brain injury, but failed in phase III clinical development for a lack of efficacy (Maas et al., 2006). This could be due, in part, to the sensitive time window between the acute insult to the brain and drug administration. Nevertheless, dexanabinol exhibits a favorable pharmacokinetic Bortezomib and dexanabinol were administered IP twice weekly (1 and 10 mg/kg). Tumor volume measurements and pharmacodynamic model fitted profiles are displayed as solid markers and dashed lines. Asterisks indicate statistical significance from control (corrected p value < 0.5), which were obtained via multiple t-tests with Holm-Sidak Bonferroni correction (GraphPad Prism v7.04).
Frontiers in Pharmacology | www.frontiersin.org January 2022 | Volume 12 | Article 817236 profile and has been shown to be safe in humans (Brewster et al., 1997;Knoller et al., 2002;Maas et al., 2006). In an effort to improve translatability over classical immortal cell lines, a combination of novel computational and experimental models was used to evaluate the potential neuroprotective effects of dexanabinol. The network-based systems pharmacology model provided an in silico projection of the neurotoxic effects of bortezomib and potential neuroprotective effects of dexanabinol on intraneuronal signaling processes. A 3D microphysiological model of peripheral nerves, although rat-derived, provides a better representation of in vivo physiology. Resident immune cells may be present in the peripheral nerve model; however, it lacks the chemotaxis and infiltration of peripheral immune cells and neuroimmune cells in the central nervous system, which have been increasingly recognized as important drivers of CIPN (Brandolini et al., 2019). Bortezomib has been shown to cause a dose-dependent decrease in neurite outgrowth in induced pluripotent stem cell (iPSC) derived neurons (Wing et al., 2017). We conducted a neurite outgrowth study using human iPSC-derived peripheral neurons (Peri.4U, Axiogenesis, Cologne, Germany); however, dexanabinol did not appear to prevent bortezomibinduced decreases in neurite outgrowth. Dexanabinol showed mixed, but promising effects in altering bortezomib-induced damage in the microphysiological nerve-on-a-chip, especially at lower concentrations (Figure 3). In addition, it restored pain thresholds that were significantly decreased in a rat model of BIPN (Figure 5), along with a partial recovery of the loss of IENF density induced by bortezomib at the end of experiment ( Figure 6). Neurophysiological studies investigate only the population of large myelinated fibers, whereas behavioral tests and IENF explore unmyelinated and small-diameter myelinated fibers (i.e., those involved in neuropathic pain). Therefore, our data indicate that dexanabinol can attenuate (although not completely prevent) the effect of bortezomib, justifying the claim that it is a possible option for bortezomib-induced neuropathic pain. These in vitro and preclinical results support the hypothesis that dexanabinol may attenuate BIPN associated pain, which was generated from a systems pharmacological analysis of peripheral neuron signal transduction.
Once the potential beneficial effect of dexanabinol was confirmed in vivo, it was critical to determine whether it interferes with the ability of bortezomib to treat myeloma. To address this concern, the interaction of dexanabinol and bortezomib was assessed in preclinical models of multiple myeloma. Dexanabinol was shown not to compromise the anticancer effects of bortezomib in these experimental models ( Figures  2B, 7C). The in vitro antagonistic effects observed in U266 multiple myeloma cells were transient and outside clinically relevant drug concentrations. Dexanabinol alone has been evaluated for the treatment of solid tumors, brain cancer, hepatocellular carcinoma, and pancreatic cancer (NCT01489826, NCT01654497, NCT02423239). Interestingly, dexanabinol appeared to enhance the anti-cancer effects of bortezomib in mice bearing MM1.S multiple myeloma xenografts ( Figure 7C). Other drugs, such as lenalidomide and dexamethasone, which exhibit anti-TNFα properties, have also been given successfully with bortezomib for the treatment of multiple myeloma.
Quantitative and systems pharmacology (QSP) has emerged as an important discipline for evaluating disease and drug effects on biological systems in a holistic manner, providing useful tools and guiding experiments throughout various stages of drug discovery and development (Sorger et al., 2011). In this study, we used a general mathematical representation of bortezomib neurotoxic effects to identify a rational polypharmacological treatment strategy for CIPN. Complex diseases with multi-factorial etiologies, such as autoimmune and neurodegenerative disease, have not been effectively managed through single-target interventions. Hence, a rational polypharmacological approach that moves beyond Paul Erhlich's magic bullet concept, a paradigm that has driven modern medicine for over a century, could lead to therapies for these unmet medical needs (Strebhardt and Ullrich, 2008;Barabasi et al., 2011). However, several study limitations must be acknowledged. First, peripheral neuropathy is a complex disorder and recapitulating the major signaling and physiological conditions are challenging. Although key intraneuronal pathways were featured in the computational model, there are other system components that may be important. As additional knowledge of the system accumulates, the signaling model could be refined to better evaluate intraneuronal toxicodynamics. Second, data availability and the translation of preclinical experiments to human therapeutic options are limited. The data used to develop and validate the systems pharmacology model were extracted from the public domain, and the original experimental conditions are often indirectly related (e.g., use of the ATF4 dataset for validation purposes). Compounds showing promise in animal experimental models of peripheral neuropathy may also fail in clinical settings. The nerve-on-a-chip model in this study utilized rat cells, and further development of microphysiological systems using human primary or iPSCderived neurons might provide additional insights into species differences in outcomes (Kramer et al., 2020). Lastly, further research is needed to test and verify the mechanisms by which dexanabinol elicits its effects on CIPN in vivo and whether these preclinical results will ultimately translate to human therapeutics. The systems pharmacology model (Supplementary Figure S1) could be expanded to incorporate additional components at the subcellular level as well as additional cell types (e.g., central and peripheral immune cells, satellite glial cells, and Schwann cells). Additionally, the model could be modified to capture quantitative relationships between system components (Eduati et al., 2017) and enable other applications, such as the ability to screen molecules with activity towards proteins in the model and support lead candidate selection. Interestingly, several components of the molecular network are also involved in central neurodegenerative diseases of the proteinopathy class, such as Alzheimer's and Parkinson's disease. For example, oxidative stress and inflammation are two hallmark characteristics of neurodegenerative disease (Fischer and Maier, 2015), and memantine, an NMDA receptor antagonist, is approved for the treatment of moderate-to-severe Alzheimer's disease. This study highlights the use of a systems pharmacology approach to reposition a drug candidate for the treatment of neuropathic pain based upon its polypharmacological profile. Owing to its favorable pharmacokinetic and safety profile in humans and preclinical efficacy, dexanabinol might represent a clinically meaningful treatment for bortezomib-induced neuropathic pain.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the co-corresponding authors.

ETHICS STATEMENT
Rat studies were reviewed and approved by the Ethics Committee of the University of Milano-Bicocca and Ministry of Health (approval numbers 0007583/19 and 842/2018-PR) and the Institutional Animal Care and Use Committee (IACUC) at Tulane University. Mouse studies were approved by the IACUC at Roswell Park Comprehensive Cancer Center.

AUTHOR CONTRIBUTIONS
All co-authors were involved in drafting the manuscript and gave approval of this version to be published. Made substantial contributions to overall study conception and design, developed and utilized systems pharmacology model, acquired in vitro cytotoxicity data, and analyzed and interpreted data, PB; Contributed to study design and acquisition of in vivo results, and analyzed and interpreted the data, CM; Acquired ex vivo nerveon-a-chip data, and analyzed and interpreted data KP; Acquired in vivo nerve conduction study data, PA and GF, and analyzed and interpreted data, PA; Assisted in obtaining in vivo model of CIPN and acquired pain data, AC, ACh, LM, EP, EB, and NO; Acquired, analyzed and interpreted IENF density data, ACh; Contributed to study design and acquisition of in vivo tumor growth dynamic data, and analyzed and interpreted data, LC; Acquired in vivo tumor growth dynamic data, and analyzed and interpreted data, WL; Acquired in vitro macrophage cytokine data, and analyzed and interpreted data, MG and TI; Contributed to study design of the in vivo tumor growth dynamic studies, and interpreted data, KL; Contributed to study design of the ex vivo nerve-on-a-chip studies, and interpreted data, MM; Made substantial contributions to overall study conception and design, and interpreted data, GC and DM.