Examining Sodium and Potassium Channel Conductances Involved in Hyperexcitability of Chemotherapy-Induced Peripheral Neuropathy: A Mathematical and Cell Culture-Based Study

Chemotherapy-induced peripheral neuropathy (CIPN) is a prevalent, painful side effect which arises due to a number of chemotherapy agents. CIPN can have a prolonged effect on quality of life. Chemotherapy treatment is often reduced or stopped altogether because of the severe pain. Currently, there are no FDA-approved treatments for CIPN partially due to its complex pathogenesis in multiple pathways involving a variety of channels, specifically, voltage-gated ion channels. One aspect of neuropathic pain in vitro is hyperexcitability in dorsal root ganglia (DRG) peripheral sensory neurons. Our study employs bifurcation theory to investigate the role of voltage-gated ion channels in inducing hyperexcitability as a consequence of spontaneous firing due to the common chemotherapy agent paclitaxel. Our mathematical investigation of a reductionist DRG neuron model comprised of sodium channel Nav1.7, sodium channel Nav1.8, delayed rectifier potassium channel, A-type transient potassium channel, and a leak channel suggests that Nav1.8 and delayed rectifier potassium channel conductances are critical for hyperexcitability of small DRG neurons. Introducing paclitaxel into the model, our bifurcation analysis predicts that hyperexcitability is highest for a medium dose of paclitaxel, which is supported by multi-electrode array (MEA) recordings. Furthermore, our findings using MEA reveal that Nav1.8 blocker A-803467 and delayed rectifier potassium enhancer L-alpha-phosphatidyl-D-myo-inositol 4,5-diphosphate, dioctanoyl (PIP2) can reduce paclitaxel-induced hyperexcitability of DRG neurons. Our approach can be readily extended and used to investigate various other contributors of hyperexcitability in CIPN.

Chemotherapy-induced peripheral neuropathy (CIPN) is a prevalent, painful side effect which arises due to a number of chemotherapy agents. CIPN can have a prolonged effect on quality of life. Chemotherapy treatment is often reduced or stopped altogether because of the severe pain. Currently, there are no FDA-approved treatments for CIPN partially due to its complex pathogenesis in multiple pathways involving a variety of channels, specifically, voltage-gated ion channels. One aspect of neuropathic pain in vitro is hyperexcitability in dorsal root ganglia (DRG) peripheral sensory neurons. Our study employs bifurcation theory to investigate the role of voltage-gated ion channels in inducing hyperexcitability as a consequence of spontaneous firing due to the common chemotherapy agent paclitaxel. Our mathematical investigation of a reductionist DRG neuron model comprised of sodium channel Na v 1.7, sodium channel Na v 1.8, delayed rectifier potassium channel, A-type transient potassium channel, and a leak channel suggests that Na v 1.8 and delayed rectifier potassium channel conductances are critical for hyperexcitability of small DRG neurons. Introducing paclitaxel into the model, our bifurcation analysis predicts that hyperexcitability is highest for a medium dose of paclitaxel, which is supported by multi-electrode array (MEA) recordings. Furthermore, our findings using MEA reveal that Na v 1.8 blocker A-803467 and delayed rectifier potassium enhancer L-alpha-phosphatidyl-D-myo-inositol 4,5-diphosphate, dioctanoyl (PIP 2 ) can reduce paclitaxel-induced hyperexcitability of DRG neurons. Our approach can be readily extended and used to investigate various other contributors of hyperexcitability in CIPN.

INTRODUCTION
Chemotherapy-induced peripheral neuropathy (CIPN) is a painful, dose-limiting side effect of chemotherapy cancer treatment that affects more than 85% of patients during treatment (Fallon, 2013) and 60% of patients 3 months post chemotherapy treatment (Seretny et al., 2014). These patients report enduring "pin and needle" paresthesias (i.e., tingling, numbness) in the peripheral nervous system (hands and feet) (Park et al., 2013). The onset of symptoms can range from 1 day to 2 years after treatment and can even persist throughout life, causing a significant decrease in the quality of life of cancer survivors (Cavaletti and Marmiroli, 2010;Speck et al., 2012). Furthermore, CIPN has been shown to result in depression, frustration, and a sense of loss of purpose among these patients (Tofthagen, 2010). To improve quality of life of cancer patients, it is imperative to find CIPN preventive agents. Since there is currently no FDA-approved treatment for CIPN, management of CIPN-induced pain includes many options such as antidepressants, anticonvulsants, anti-inflammatory, and opioid therapies.
Neuropathic pain can be characterized by changes in firing patterns, calcium signaling, mitochondrial dysfunction, axonal transport dysfunction, and inflammation, among others (Meacham et al., 2017). Increases in spontaneous firing has been strongly associated with allodynia but hyperalgesia may be related to reduced threshold of triggered firing or enhanced firing of putative small pain sensing DRG neurons (Devor and Seltzer, 1999). Hyperexcitability from increased spontaneous firing involves voltage-gated sodium channels Na v 1.6-1.9 (Lou et al., 2005) and potassium channels (Du and Gamper, 2013), among others. Although there is some role of the central nervous system (CNS) in CIPN, treatments targeting CNS-based pain pathways have not been sufficient to reduce CIPN (Chu et al., 2015;Fukuda et al., 2017). Mainly, CIPN has been studied in the dorsal root ganglia (DRG) sensory neurons in the peripheral nervous system (PNS). DRG neurons in the PNS relay sensory signals to neurons in the spinal cord and eventually to the CNS, thus allowing the different subpopulations of DRG neurons to respond to different nociceptive stimuli including mechanical, thermal, and chemical (Dubin and Patapoutian, 2010). DRG neurons are more susceptible to chemotherapy agents than the CNS neurons because DRG neurons do not have an extensive neurovascular barrier to limit drug entry (Lees et al., 2017;Livni et al., 2019). Thus, we chose DRG neurons to be our model system for this study.
A myriad of alterations in various pathways have been linked to CIPN including voltage-gated ion channels, in calcium signaling, in fast axonal transport, and in occurrence of oxidative stress and inflammation (Carozzi et al., 2015). Several chemotherapy agents, such as vincristine, paclitaxel, and oxaliplatin have been suggested to cause CIPN (Carozzi et al., 2015). In this study, we focus on the chemotherapy agent paclitaxel (brand name Taxol). Paclitaxel is a microtubulebinding cancer agent used for several solid tumor cancers such as breast, ovarian, and lung. The conjecture concerning the paclitaxel-induced CIPN (PIPN) mechanism is that it reduces axonal transport of mRNA which can lead to axonal degeneration, alters expression of membrane ion channels, and induces inflammation and oxidative stress (Carozzi et al., 2015). Several agents have been tested in clinical trials, but their ability to prevent PIPN is still unclear (Argyriou et al., 2008). Since it can be difficult to control and examine multiple events in an experimental setting, a mathematical model can be employed to investigate mechanisms by exploring many channels in a controlled manner to observe how the system behavior changes upon any external influences. In this study, we analyze the role of voltage-gated sodium and potassium ion channels in paclitaxel-induced hyperexcitability using mathematical modeling and bifurcation theory, then support our results using in vitro recordings. Specifically, our mathematical model is representative of a simplified small DRG neuron model consisting of two sodium channels (Na v 1.7 and Na v 1.8), two potassium channels [delayed rectifier (KDR) and Atype transient (KA)], and one leak channel. We included these channels because they were deemed to be most prominent in small DRG neurons (Gold et al., 1996;Sheets et al., 2007;Choi and Waxman, 2011). Here, we report the development of a mathematical modeling-based approach, which is then supported by experimental findings in cell culture-based model of CIPN.

Mathematical Model Description and Equations
The mathematical model here is representative of a small DRG neuron which is simplified to be a single compartment cylindrical model. It consists of sodium channels Na v 1.7 and Na v 1.8, delayed rectifier (KDR) and A-type transient (KA) potassium channels, and a leak channel. Equations and parameter values used in this model were from Choi and Waxman (2011). The main equation for membrane voltage is written as: where, I ext is the external applied current, i 1.7 , i 1.8 , i KDR , i KA , i l are specific ionic currents due to Na v 1.7, Na v 1.8, KDR, KA, and leak channels. Model parameters C is the specific membrane capacitance and A is the area, and the variables V is the membrane voltage and t is time. These ionic currents are written as following: where, g i (i=1.7, 1.8, KDR, KA, l) are maximal conductances and are constants. Parameters E Na , E K , and E l are the equilibrium ion potentials. All the activation and inactivation gating variables x  (Hodgkin and Huxley, 1952): The expressions of x ∞ and τ x , and the parameter values are specified below. All the equation forms have been extracted from literature (Choi and Waxman, 2011). The model parameter values and the corresponding references are mentioned in Table 1. Below, we describe the equations for each of the voltagegated ion channels. The equations for Na v 1.7 are written as the following: For x = m 1.7 , h 1.7 , s 1.7 , , and, The equations for Na v 1.7 were extracted from Sheets et al. (2007) and Choi and Waxman (2011). Variables m 1.7 is the activation gating variable, h 1.7 the fast-inactivation gating variable, and s 1.7 the slow-inactivation gating variable. The equations for Na v 1.8 are written as the following: , and, The equations for Na v 1.8 were obtained from Sheets et al. (2007) and Choi and Waxman (2011). Variables m 1.8 is the activation gating variable and h 1.8 is the inactivation gating variable. The equations for KDR are written as the following: where, α n KDR = 0.001265 × 10 if V = − 14.273, , and, The equations of KDR channel were obtained from Schild et al. (1994). Variable n KDR is the activation gating variable. Inactivation gating variables are not present.
The equations for KA are written as the following: , and, The equations for KA channel were obtained from Sheets et al. (2007). Variables n KA is the activation gating variable whereas h KA is the inactivation gating variable. The timescales of all these activation and inactivation variables for each of the aforementioned channels can be found in Verma et al. (2020c).

Primary Cell Culture
Dorsal root ganglia (DRG) neurons were extracted from wildtype Sprague-Dawley rat pups 7-14 days old. Age of rats were chosen based on a previous CIPN study (Li et al., 2018) since cultured postnatal DRG neurons share many characteristics of mature DRG neurons (Melli and Hake, 2009). Animals were maintained in the norovirus-negative facility of the Centrally Managed Animal Facilities at Purdue University. They were housed in at a constant temperature and humidity on a 12:12 light-dark cycle (lights on 06:00-18:00) with ad lib access to food and water according to as approved by the Purdue Animal Care and Use Committee and the Institutional Animal Care and Use Committee (IACUC) and was conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Pups were placed on a paper towel on ice for 2 min before decapitation. The spinal cord was extracted and DRG neurons were dissected and placed in complete saline solution (CSS). Cells were prepared by centrifuging at 2.5×g for 30 s then adding 1.5 mg/mL of collagenase A in CSS with 0.05 mM EDTA. After rotating in the 37 • C incubator for 20 min, cells were spun at 2.5×g for 30 s. 1.5 mg/mL collagenase D and 30U papain in CSS were added and then placed in the incubator rotator for 20 min. Cells were spun at 2.5×g for 3 min. Primary DRG neurons were triturated in 1 mL of 0.15% trypsin inhibitor and 0.15% bovine serum albumin (BSA) in Dulbecco's Modification of Eagle's Medium (DMEM) medium with 10% fetal bovine serum (FBS) (Corning, Corning, NY, USA) and then spun again at 2.5×g for 3 min then filtered through a 40µM filter before seeding on a plate coated with laminin and poly-d-lysine. DRG neurons from one pup were divided into six wells. Culture media was changed after 2 days to NbActiv4 recording media (BrainBits, Springfield, IL, USA). Four days after seeding, culture was recorded.

Micro/Multielectrode Array (MEA)
Firing properties were recorded using a Maestro Pro (Axion BioSystems, Atlanta, GA, USA). Twelve well plates of 64 electrodes per well were used for culture. Plates were coated the day of seeding by incubating with poly-d-lysine for 2 h, washing with sterile milli-q water three times, then incubating with laminin for 1 h. All recordings were at 37 • C. The plate was recorded before treatment and 24 h after treatment. Then 200 µL of 1.0 µM capsaicin was added and recorded for 2 min as a positive control. Analysis was performed using the manufacturer's software, Axion BioSystems Integrated Studio (AxIS) and NeuroMetric Tool. An electrode (n = 1) was considered active if there were more than two action potentials in the baseline, response to buffer, or response to capsaicin. Mean firing rate (Hz) was calculated for active electrodes. Chili pepper compound capsaicin was added into the system to trigger a sensory response. Fold change was calculated from natural logarithm of firing rate in line with literature (Atmaramani et al., 2020;Negri et al., 2020). Percent change described in Table 2 is the number of spontaneous firing neurons after treatment minus the number before treatment over number before treatment. Recordings were compiled from different cultures extracted from different animals on different days. Although one electrode may detect more than one neuron, this situation is unlikely to occur to a notable level, due to low basal firing and low plating density of primary DRG sensory neurons in the MEA experiment (Yang et al., 2016(Yang et al., , 2018.

Statistical Analyses
A Shapiro-Wilk test determined that the data was not normallydistributed so Krushal-Wallis and Mann-Whitney U-tests were used to determine statistical significance of treatment differences. MEA results are presented as mean (median, range, sample size, p-value). One-sided chi-square test was used to analyze the increase in firing rate from baseline to treatment. Also, this test was used to analyze the change in qualitative behavior of neurons by comparing the number of neurons with a firing rate of zero in the baseline and with a > 0 firing rate after treatment. P-values were reported at *p ≤ 0.05, **p ≤ 0.01, and ***p≤0.001. GraphPad Prism 8.3.0 was used to determine statistical significance.

Model Simulations
One indicator of peripheral neuropathy is spontaneous firing which implies repetitive firing of action potentials for I ext = 0. In this work, we explore the parameters that can lead to spontaneous firing and can be potentially impacted by paclitaxel. We vary the maximal ion conductances to explore whether they can induce spontaneous firing since current literature indicates that paclitaxel can manipulate the expression of voltage-gated ion channels (Zhang and Dougherty, 2014). A bifurcation analysis of this model with I ext as the bifurcation parameter can be found elsewhere (Verma et al., 2020c). First, we performed dynamic simulations for different parameter values. The initial conditions correspond to the stable steady state solution obtained for the parameter values mentioned in Table 1. Figure 1 demonstrates how the voltage dynamics vary upon increasing g 1.8 (g 1.8 was chosen for illustration purposes), without and with noise (following Rho and Prescott, 2012 for including noise). Without noise, for low values of g 1.8 , the system settles down to a steady state, shown in Figure 1A. For higher values of g 1.8 , mixed-mode oscillations (MMOs) are observed, shown in Figure 1B. MMOs consist of both small amplitude (subthreshold) and large amplitude (action potential) oscillations. For higher values, continuous firing of action potentials is observed, shown in Figure 1C. In the next section, we focus on the switch from steady state to continuous firing by treating different channel conductances as bifurcation parameters. We have majorly focused on the model without noise and have only demonstrated how simulations vary upon introducing noise. Mixed-mode oscillations and the impact of noise were investigated in detail in another work (Verma et al., 2020b).

Bifurcation Analysis
XPPAUT is used to perform preliminary bifurcation analysis, and confirmation of results is done using MATCONT. Parameter settings for XPPAUT and MATCONT are mentioned in section 2. As mentioned before, bifurcation diagrams are generated by setting the maximal conductance as the bifurcation parameter, one by one, for each voltage-gated ion channel involved in this model (Na v 1.7, Na v 1.7, KDR, and KA).

One-Parameter Continuation
Firstly, we performed one-parameter continuations to find bifurcation points which could separate steady state from MMOs, and MMOs from continuous firing of action potentials. A partial bifurcation diagram is shown in Figure 2. As seen in Figures 2A,D, no bifurcation points are generated upon varying g 1.7 and g KA . A single red line, representing stable steady state solutions, is observed. On the contrary, bifurcation points are observed in Figures 2B,C. In Figure 2B, a subcritical Hopf bifurcation point (HB) is detected upon increasing g 1.8 . Beyond this point, steady state solutions become unstable, shown by the black branch. Two turning points/limit points (LP 1 , LP 2 ) are also detected on the black branch of unstable steady state solutions. Since the Hopf bifurcation point is subcritical, unstable periodic solutions emanate from it, as shown by the blue branch. This branch first turns at a cyclic limit point (CLP 1 ) and then meets the unstable steady state branch, indicating a homoclinic orbit. This turning at CLP 1 is not obvious from the figure, however, it can be observed upon zooming into the branch. Upon moving in the backward direction starting from a large value of g 1.8 resulting in stable periodic solutions, a stable periodic solution branch is generated, indicated in green, finally leading to a cyclic limit point (CLP 2 ) beyond which the periodic solutions become unstable, indicated in blue. This unstable periodic branch abruptly ends due to the period of the branch increasing substantially, indicating that it may tend toward a period-infinity solution, as shown in Figure 2E. The stable periodic solution branch indicates the spontaneous firing parameter regime. The frequency of spontaneous firing increase with g 1.8 , as shown in Figure 2E. MMOs are found in some region between the Hopf bifurcation point (HB) and the cyclic limit point (CLP 2 ), shown by the shaded pink region. Moreover, as was observed in Rho and Prescott (2012), we noted that the parameter range of MMOs will become wider upon introducing noise in the model. We have not explored the bifurcations in the MMOs regime; however, a detailed discussion for a similar situation is given in Verma et al. (2020b).
A similar, although horizontally flipped, bifurcation diagram is generated with g KDR as the bifurcation parameter, shown in Figure 2C. Upon decreasing g KDR , a subcritical Hopf bifurcation point (HB) is detected with unstable periodic solutions emanating from it. These unstable periodic solutions also intersect the unstable steady state branch, indicating a homoclinic orbit. Similar to Figure 2B, a stable periodic solution branch (indicated by green circles) is detected as well which becomes unstable after a cyclic limit point (CLP 2 ). As shown in Figure 2F, the unstable periodic solution again seems to tend toward a period-infinity solution. Moreover, the frequency of spontaneous firing decreases with increase in g KDR .
These bifurcation diagrams indicate that manipulating g 1.8 or g KDR can induce spontaneous firing, while manipulating g 1.7 or g KA will not decrease hyperexcitability. Therefore, to reverse hyperexcitability, Na v 1.8 and KDR channels should be targeted. We targeted these channels in the case of paclitaxel-induced hyperexcitability, in a primary rodent DRG neuron culture, the results of which are described in the section 3.4.

Two-Parameter Continuation
We also performed two-parameter continuation in order to explore the combinational effects of these conductances. To this end, we observed changes in the detected bifurcation points upon altering another maximal conductance. In particular, we performed continuation of HB and CLP 2 . These results are shown in Figure 3.
As seen in Figures 3A,B, g 1.7 and g KA do not substantially impact the bifurcation points even in combination with g 1.8 . Decreasing g KA can shift the cyclic limit point of g 1.8 to the right, as seen in Figure 3B, which implies that the MMOs regime will become wider. In both the cases, the bifurcation points vary within a narrow range of g 1.8 . Similarly, Figures 3C,D show that g 1.7 and g KA do not impact the bifurcation points substantially even in combination with g KDR . In these two cases, an increase in g 1.7 and g KA can lead to a region of bistability where stable steady state and continuous firing of action potentials solutions coexist. However, the bifurcation points only vary within a narrow range of g KDR . Upon varying g 1.8 and g KDR together, the bifurcation points vary linearly, as shown in Figure 3E. This indicates that decreasing g 1.8 and increasing g KDR can eliminate spontaneous firing.

Effect of Paclitaxel
Current literature suggests that paclitaxel can impact gene expression of various voltage-gated ion channels (Zhang and Dougherty, 2014;Aromolaran and Goldstein, 2017). However, it is not known whether the impact is direct or indirect. Evidence indicates that paclitaxel impacts inflammatory cytokines, which can subsequently manipulate the ion channels (Aromolaran and Goldstein, 2017). For example, these inflammatory signals can increase sodium current (Wang et al., 2012). Moreover, a sigmoidal shaped dose-dependent relation is observed between paclitaxel and macrophage IL-12, as seen in Figure 3 in Mullins et al. (1999). Based on this evidence, we assumed a Hill's kinetics type relation between paclitaxel and ion channel maximal conductances. Hill's kinetics are widely used to model doseresponse curves. We assume that the conductances will vary as a function of paclitaxel dosage. Moreover, we assume that paclitaxel will lead to an increase in maximal conductance of both the sodium channels, while it would lead to a decrease in maximal conductance of both the potassium channels since these cases would lead to spontaneous firing. We assume that all the conductances are impacted even though it may not lead to a change in spontaneous firing, as shown in the previous section. Therefore, we consider the following relationships: g 1.7,new = g 1.7 + (G Na,max − g 1.7 ) [P] h n where, h n is Hill's coefficient, [P] is the paclitaxel concentration (in nM), k 0.5 is the half maximal effective concentration, g i (i=1.7, 1.8, KDR, KA) is the original maximal conductance value, g i,new is the updated maximal conductance value from the above equation. The upper and lower limit of maximal conductances is represented by G Na,max and G K,min , respectively. Depending on G j (j=Na, max, K, min) and h n , the relationship between [P] and g i,new (i=1.7, 1.8, KDR, KA) will vary, as shown in Supplementary Figures 1A-D. Increasing or decreasing G j increases or decreases the maximal conductance parameter (g i,new ) range covered upon varying [P], respectively. Large h n creates sigmoidal curves. Decreasing h n makes the curve seem exponential. In the remaining text, representative values for G j and h n were taken as listed in Table 1 and the relationship is demonstrated in Figure 4D. In Supplementary Figure 2, a sensitivity analysis with respect to these parameters was conducted, which shows that the qualitative behavior does not vary much upon varying these parameter values.
Following this, we performed numerical bifurcation analysis with paclitaxel concentration [P] as the bifurcation parameter. A partial bifurcation diagram is shown in Figure 4A which reveals that upon stable steady state solution continuation, a subcritical Hopf bifurcation point (HB 1 ) is found beyond which the solutions become unstable. Unstable periodic solutions emanate from this Hopf bifurcation point. Upon continuation of the unstable steady state solution branch in black, a supercritical FIGURE 3 | Two parameter continuations performed for the Hopf bifurcation point HB and cyclic limit point CLP 2 . (A) Continuation plot for g 1.7 vs. g 1.8 show that the bifurcation points generated by keeping g 1.8 as the bifurcation parameter do not shift upon varying g 1.7 . There is a range of g 1.8 where MMOs can be observed. (B) The HB of g 1.8 bifurcation diagram do not shift upon varying g KA . CLP 2 shift rightwards upon decreasing g KA . This implies that the MMOs region will be wider in this case. (C) Bifurcation points of g KDR do not shift upon varying g 1.7 . For a narrow range of g KDR , MMOs and bistability between stable steady state and continuous firing are observed. (D) The HB of g KDR do not shift upon varying g KA . CLP 2 shifts leftwards upon decreasing g KA . This implies that the MMOs region will become narrower in this case. A region of bistability is also observed for a narrow range of g KDR . (E) A linear combinational effect is seen between g 1.8 and g KDR . Note that the thin gap between stable steady state and continuous firing regimes is the MMOs region.
Hopf bifurcation point (HB 2 ) is found, beyond which the steady state solutions become stable again. Stable periodic solutions (indicated in green) emanate from this point. The paclitaxelinterval between these two Hopf bifurcation points constitutes the spontaneous firing regime. The subinterval where the periodic and the steady state branch are unstable (between PD and CLP 4 ) corresponds to the MMOs regime. MMOs of this model have been studied in some detail in Verma et al. (2020b); a detailed investigation of the MMOs shown in Figure 4A was beyond the scope of this work. Stable periodic solutions with a small amplitude arise from the supercritical Hopf bifurcation point which turn unstable after a cyclic limit point (CLP 4 ). The unstable periodic solutions finally become stable after a subcritical period doubling bifurcation point (PD) when going in the direction of decreasing [P]. The stable periodic solutions become unstable again after a cyclic limit point (CLP 2 ).
The frequency of firing in the first stable periodic solutions regime (between CLP 2 and PD) is shown in Figure 4C. It is shown that upon increasing paclitaxel concentration, frequency of firing first increases, and then decreases after reaching a maximum firing rate. Beyond the PD point, the frequency of firing decreases further since the solutions are of MMOs type.

Experimental Validation Results
A high-throughput way to continuously record neuronal firing patterns is to use multielectrode array (also referred to as microelectrode array, MEA). This system is capable of recording of extracellular voltage potentials with millisecond temporal resolution of neurons in cultures grown on a 768 array of electrodes up to 96 well format, which makes it high-throughput.
To support the findings of the mathematical modeling of the effect of paclitaxel dose on hyperexcitability, we measured the firing rate for different doses of paclitaxel. Low doses (10 nM) and high doses (1 µM) of 24-h paclitaxel administration caused lower firing rate than 250 nM paclitaxel as expected (see Supplementary Statistics File). Thus, we decided to use 250 nM for the dose of the subsequent experiments. Na v 1.8 blocker A-803467 and KDR enhancer Lalpha-phosphatidyl-D-myo-inositol 4,5-diphosphate, dioctanoyl (PIP 2 ) when administered separately together with paclitaxel reduces the number of spontaneous firing neurons ( Table 2) and the firing rate ( Figure 5A). Similarly, the representative heat maps of firing frequency reveals the same trend qualitatively ( Figure 5C). Specifically, the natural log of firing rate fold change (treatment/baseline) increases with paclitaxel addition to a mean of 0.15 (median = 0.016, range = −2.97 to 3.26, n = 442, p < 0.0001 compared to media). Na v 1.8 blocker A-803467 decreases the natural log of firing rate to a mean of 0.08 (median = 0.001, range = −2.40 to 2.50, n = 338, p = 0.0451 compared to paclitaxel). Similarly, KDR enhancer PIP 2 reduced paclitaxelinduced hyperexcitability to 0.08 (median = 0.00, range = −3.34 to 2.86, n = 266, p < 0.0001 compared to paclitaxel).
A chi-square analysis revealed that paclitaxel-treated wells had more neurons whose firing rate increased from baseline compared to the media control (χ 2 = 22.47, z = 4.74, p < 0.0001). Wells treated with Na v 1.8 blocker A-803467 had less neurons whose firing rate increased from baseline as compared to those with only paclitaxel treatment (χ 2 = 5.85, z = 2.42, p = 0.0078). Similarly, wells treated with KDR enhancer PIP 2 had less neurons whose firing rate increased from baseline compared to those with only paclitaxel treatment (χ 2 = 86.55, z = 9.30, p < 0.0001). We also compared the MEA results with our bifurcation results from Figures 2, 4 which showed that increasing g 1.8 or [P] to some extent, or decreasing g KDR can lead to spontaneous firing. To this end, we analyzed the number of neurons that did not fire in the baseline, but fired in the treatment, which we considered FIGURE 5 | Analysis of multielectrode array (MEA) recording summary shows amelioration of hyperexcitability after treatment of A-803467 (Na v 1.8 blocker) and PIP 2 (KDR enhancer). (A) Firing rate reveals a significant increase in paclitaxel firing from media control (p < 0.0001), a significant decrease from paclitaxel when A-803467 or PIP 2 are administered concurrently with paclitaxel (p = 0.0451 and p < 0.0001, respectively) (n = 275 for media, 442 for paclitaxel, 338 for paclitaxel and A-803467, 266 for paclitaxel and PIP 2 ). Box and whisker plots include a box at the 25 th and 75 th percentiles, whiskers extend to the 10 th and 90 th percentiles, and the line is at the median. (B) Dynamic neurons are described as neurons that had a firing rate of zero in the baseline with an increased firing rate after treatment. Bottom bar (black or colored) is the percent of neurons that had a firing rate of zero in baseline with an increased firing rate after treatment. Top bar (gray) is the percent neurons that did not fall in this category. (C) Heatmap of representative MEA recordings with firing frequency of each active electrode color-coded with warm colors (red, orange, yellow) representing high firing frequency (white = 10 Hz) and cool colors (green, blue) representing low firing frequency (black = 0 Hz). Each circle represents a spontaneously firing neuron within the 8 × 8 electrode array. Top row is baseline at time zero before treatment is added. Bottom row is 24 h after treatment was added. Asterisks denote statistical significance from Mann-Whitney U-test (A) or chi-square (B) (*P < 0.05, ***P < 0.001).
as a change in the qualitative behavior of the dynamics of these neurons ( Figure 5B). We refer to these neurons as dynamic neurons. Paclitaxel-treated wells had an increased number of dynamic neurons compared to the media control (χ 2 = 10.28, z = 3.21, p = 0.0007). The Na v 1.8 blocker A-803467 treated group had a decrease in dynamic neurons compared to paclitaxel treated group (χ 2 = 3.46, z = 1.86, p = 0.03) as did the KDR enhancer PIP 2 treated group (χ 2 = 31.77, z = 5.64, p < 0.0001). These results support the hypothesis that blocking Na v 1.8 or enhancing KDR can reduce paclitaxel-induced hyperexcitability.

DISCUSSION
CIPN is a debilitating experience for cancer patients with no current established methods of preventing or treating it due to minimal understanding of its pathophysiology (Park et al., 2013). In this paper, we apply a mathematical approach using numerical bifurcation theory to understand the role of sodium channel Na v 1.7, sodium channel Na v 1.8, delayed rectifier potassium channel, and A-type transient potassium channel in CIPN in a putative small DRG neuron model. Maximal conductances were kept as bifurcation parameters to identify which of these channels can induce spontaneous firing (one of many indicators of peripheral neuropathy). We further used MEA experiments to support our findings.
Using bifurcation theory, we found that, increasing g 1.8 and decreasing g KDR can induce spontaneous firing (see Figure 2) in a simplified DRG neuron model. The effect may be aggravated in combination, as seen from the two parameter plot in Figure 3E. Our results indicate that a Na v 1.8 blocker should reduce spontaneous firing which supports the role of Na v 1.8 in contributing to the increased excitability in peripheral neuropathy (Xiao et al., 2019;Zhang et al., 2019). The significance of blocking Na v 1.8 in PIPN was also observed in our MEA experiment in which Na v 1.8 blocker A-803467 had a neuroprotective effect on paclitaxel-induced increased firing rate when administered at the same time as paclitaxel in primary DRG neuron culture (Figure 5). Similarly, KDR was indicated in our model and MEA findings to be involved with hyperexcitability which is supported by literature (Du and Gamper, 2013). Although these results are for specific dosages of A-803467 and PIP 2 to observe the change in excitability, this electrophysiology data supports the trends found in the mathematical model analysis. Depending on the amount of A-803467 and PIP 2 , hyperexcitability should change. Interestingly, Na v 1.8 and KDR channels were found to be sensitive even with I ext as the primary bifurcation parameter (Verma et al., 2020b). It will be of interest to investigate their protective effects on CIPN due to other chemotherapy agents such as vincristine and oxaliplatin, as well.
We assume a Hill's type kinetics for the effect of paclitaxel on the ion channels. Based on this, a partial bifurcation diagram is generated, treating paclitaxel concentration as the bifurcation parameter. This bifurcation diagram indicates that spontaneous firing should arise for a mid-range of paclitaxel dosage. Moreover, firing rate should first increase and then decrease, seen from the frequency diagram in Figure 4C. A similar trend is seen in the MEA recordings, shown in the Supplementary Statistics File. Firing rate is larger for the middle value of paclitaxel concentration in the range considered here. The frequency diagram ( Figure 4C) may not reflect the actual reason for this trend. It may be that the cells died at a higher concentration, or they may not fire because of manipulation of the ion channels as shown by a mathematical relationship with paclitaxel. Further investigation is required to establish the reason behind this observation with certainty.
The parameters for the relationship between paclitaxel and ion channel maximal conductances are also assumed (see Table 1). In the Supplementary Figure 2, we have shown how the behavior of the model varies upon changing these parameters. Upon varying h n , G Na,max , or G K,min , the qualitative behavior of the diagram does not change for the range of parameters considered here; stable steady states, MMOs, and continuous firing regimes are observed. Assuming that Hill's kinetics is a reasonable relationship between paclitaxel and ion channel conductance, defined upper and lower limits of how much the conductances can vary may exist. These parameters can be estimated by recording currents due to each of these ion channels for different doses of paclitaxel, using patch-clamping experiments. Lastly, the value of k 0.5 we used in our mathematical model (500 nM) is different from the paclitaxel dose amount used for the MEA experiment (250 nM). This is because we estimated the value from literature which was based on the IC50 value for iPSCderived human neurons and clinical data based on paclitaxel's toxicity on cancer cells (Rowinsky et al., 1999;Rana et al., 2017). However, as seen in Figure 4C, maximum spontaneous firing is observed around 300 nM. Since we only intend to perform a qualitative comparison between the bifurcation diagram ( Figure 4A) and MEA paclitaxel dose trend (included in Supplementary Statistics File), the value of k 0.5 assumed is not important. A different value of k 0.5 will shift the complete bifurcation diagram, however, the qualitative structure of the diagram will remain the same.
The mathematical model that we considered is a minimal model representing dynamics of one type of peripheral neurons, putative small pain-sensing DRG neurons. This model provides a framework for future studies of other channels, pumps, and exchangers in hyperexcitability due to CIPN. In vivo, many more ion channels are present in these DRG neurons and they can be added to the model in future studies. For example, transient receptor potential channels (TRP) have also been suggested to be involved in CIPN (Hara et al., 2013;Chukyo et al., 2018). In addition, calcium channels can play a role in CIPN (Schmitt et al., 2018). These molecules will be excellent candidates for future studies, which may provide additional insights to CIPN. More detailed mathematical models have been developed previously and can be used for this purpose (Mandge and Manchanda, 2018). Paclitaxel can also induce cytosolic calcium oscillations (Boehmerle et al., 2006), which can again be analyzed using bifurcation theory. Another factor to add in future models is the impairment of axonal transport that occurs due to paclitaxel and its effect on microtubules (Nicolini et al., 2015), which can be modeled using cable equations (Holmes, 2013). In addition, a multicompartment modeling approach can be undertaken to investigate how spike propagation across the T-junction of the DRG neuron is regulated by the ion channels (Gemes et al., 2013;Du et al., 2014;Sundt et al., 2015). As mentioned previously, the effect of paclitaxel on ion channels may be indirect and due to other inflammatory cytokines (Aromolaran and Goldstein, 2017). In the future, this can also be included in the next phase of modeling. Lastly, our model consists of only small DRG neurons. However, paclitaxel impacts medium and large DRG neurons as well (Zhang and Dougherty, 2014). It will be of interest to evaluate models representing different subtypes of DRG neurons and investigate the role of differentially expressed channels specific to each subtype. The MEA experiment recorded firing of all DRG neuron subtypes. The blocker and enhancer used in the MEA experiments may have regulated ion channels besides those in small DRG neurons. Similarly, paclitaxel may have impacted the firing rate of different DRG neuron subtypes as well as other supporting cells such as glia. This will require further investigation in the future by recording DRG neuron subtype-specific firing rate. While additional compartments will be needed to make the model more comprehensive, its interpretation becomes challenging with increasing complexity. With the current model, it is difficult to explore how the MMOs arise due to the multiple time scales. A model dimension reduction method such as that mentioned in Rubin and Wechselberger (2007) will be required for an analytical investigation.
DRG neuron firing has been correlated with neuropathic pain in vitro (Yang et al., 2016(Yang et al., , 2018, thus analyzing changes in firing patterns allows for a surrogate to study "pain in a dish." Recent advances in neuroelectrophysiology technology have allowed for more temporal and spatial dynamic data collection on live cells. We chose a non-invasive method to measure the electrophysiological properties (spontaneous/non-evoked firing rate) of a network of neurons through microelectrode array (MEA) recording. Therefore, we used MEA to support our model outputs by recording the neuronal activity of the effect of paclitaxel on cultured DRG neurons, as shown in Figure 5. Although MEA cannot provide specific electrophysiological firing patterns of non-spontaneously firing neurons to detect bifurcation points ex vivo, this method was chosen because it is non-invasive, high-throughput, and can access neurons in their regular culture medium which can be more physiologically relevant. Even though a change in firing rate does not demonstrate that the system undergoes a bifurcation, one can use the change in the number of spontaneously firing neurons and a change in firing rate as an indicator of a qualitative change in the dynamics of the neurons. Using MEA, we can investigate how this changes upon introducing paclitaxel and the modulators. As seen in Table 2, the number of spontaneously firing neurons increases upon introducing paclitaxel and decreases upon introducing the modulators, indicating a qualitative change in the dynamics of an increasing number of neurons. Moreover, the number of neurons that spontaneously fired after treatment was significantly higher with paclitaxel treatment alone as compared to paclitaxel treatment together with modulators ( Figure 5B). In this way, MEA can be indirectly used to support our bifurcation analysis findings. In the future, it will be of interest to explore a strict evidence for the bifurcations by recording dynamics of a single neuron. Moreover, it will also be of interest to match the spiking patterns from the model, in particular the mixed-mode oscillations, to the recordings from a single neuron.
Our bifurcation analysis identified that the maximal conductances of Na v 1.8 and KDR are involved with the induction of spontaneous firing, associated with paclitaxel-induced hyperexcitability. The results from this model demonstrated that paclitaxel affects hyperexcitability in a concentration-dependent manner, and that Na v 1.8 and KDR conductances can impact this hyperexcitability. We propose that one of the mechanisms behind PIPN is paclitaxel regulating expression of these ion channels-which are essential for hyperexcitability to occur-in a concentration-dependent matter, either directly or indirectly. This approach, along with patient-specific pharmacokinetics of paclitaxel, can be applied to the clinic by obtaining the patient's electrophysiological profile using patch clamp or MEA and by using those values as parameters in our model to determine the optimal dose of paclitaxel based on the individual, and for examining acute vs. chronic onset of PIPN. Moreover, treatments can be designed to specifically block Na v 1.8 or enhance KDR conductance to reduce hyperexcitability caused by paclitaxel since we have identified the corresponding channels as key contributors to the hyperexcitability related to PIPN.

DATA AVAILABILITY STATEMENT
The XPP code for this model (with and without paclitaxel) can be found in the ModelDB database (http://modeldb.yale.edu/ 264591).

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional (Purdue) Animal Care and Use Committee.

AUTHOR CONTRIBUTIONS
PV, AK, DF, and DR designed and supervised the numerical bifurcation study. PV generated the bifurcation results. ME and YY designed the experimental study and ME performed the MEA experiments. PV and ME wrote the manuscript. Everyone participated in designing this study and reviewing the manuscript.