Computational analysis of a 9D model for a small DRG neuron
Parul Verma, Achim Kienle, Dietrich Flockerzi, Doraiswami Ramkrishna
CComputational analysis of a 9D model for a small DRG neuron
Parul Verma Achim Kienle , Dietrich Flockerzi , ∗ Doraiswami Ramkrishna Davidson School of Chemical Engineering, Purdue University, USA Max Planck Institute for Dynamics of Complex Technical Systems, Magdeburg, Germany Otto von Guericke University, Magdeburg, Germany ∗ [email protected] Abstract
Small dorsal root ganglion (DRG) neurons are primary nociceptors which are responsible for sensing pain.Elucidation of their dynamics is essential for understanding and controlling pain. To this end, we presenta numerical bifurcation analysis of a small DRG neuron model in this paper. The model is of Hodgkin-Huxley type and has 9 state variables. It consists of a Na v v v g . ). Numerical bifurcation analysis shows that depending on g . and the externalcurrent, different parameter regions can be identified with stable steady states, periodic firing of actionpotentials, mixed-mode oscillations (MMOs), and bistability between stable steady states and stable periodicfiring of action potentials. We illustrate and discuss the transitions between these different regimes. Wefurther analyze the behavior of MMOs. As the external current is decreased, we find that MMOs appearafter a cyclic limit point. Within this region, bifurcation analysis shows a sequence of isolated periodicsolution branches with one large action potential and a number of small amplitude peaks per period. Fordecreasing external current, the number of small amplitude peaks is increasing and the distance between thelarge amplitude action potentials is growing, finally tending to infinity and thereby leading to a stable steadystate. A closer inspection reveals more complex concatenated MMOs in between these periodic MMOsbranches, forming Farey sequences. Lastly, we also find small solution windows with aperiodic oscillations,which seem to be chaotic. The dynamical patterns found here as a function of different parameters containinformation of translational importance as their relation to pain sensation and its intensity is a potentialsource of insight into controlling pain. Neurons display a variety of rich dynamics such as repetitive firing of action potentials, bursting, mixed-mode oscillations, and bistability. The diversity of dynamics displayed by a multitude of neurons has ledseveral researchers to bifurcation theory to understand the transition from one dynamical pattern to the otherfor more than 40 years [1–4]. Starting with the analysis of low dimensional models such as Hodgkin-Huxleyand Fitzhugh-Nagumo equations, it has recently been used for higher dimensional models such as a 14D1 a r X i v : . [ q - b i o . N C ] J a n odel of a pyramidal cell [5], as well. In this paper, we employ numerical bifurcation analysis to understandthe dynamics of a 9D model of a small dorsal root ganglion (DRG) neuron.Small DRG neurons are primary nociceptors, i.e., they are responsible for sensing pain [6]. From atheoretical point of view, pain corresponds to repetitive firing of action potentials [7, 8]. To understandhow pain can be controlled, it is therefore essential to determine how the transition to periodic firing ofaction potentials depends on the physiological parameters and how these parameters can be manipulated ina suitable way.While limited numerical [9–11] and extensive experimental [12–15] studies have been executed for thistype of neuron, a detailed bifurcation analysis has not been undertaken so far. The importance of usingbifurcation theory to understand pain is emphasized in the works of [16] and [17] where 2D and 3D modelsof an afferent sensory neuron were analyzed with regard to neuropathic pain and, subsequently, bifurcationtheory aided in finding parametric regions of pain and no-pain. Previous work on a model of a small DRGneuron [18] also illustrates the utility of bifurcation theory for understanding pain. In that paper, geneticmutations in sodium channels that are associated with pain sensation were also investigated.In the present paper, we use the aforementioned theory extensively to find the bifurcations explainingthe excitability patterns of this model. We perform both one-parameter and two-parameter continuation ofmodel solutions, with external current as the primary bifurcation parameter and maximal conductance ofone of the voltage-gated ion channels as the secondary parameter. Here, we find different solution regimesconsisting of stable steady states, periodic firing of action potentials, and mixed-mode oscillations (MMOs).The latter are periodic or aperiodic oscillatory solutions consisting of small amplitude (subthreshold) andlarge amplitude (action potential) peaks. They have been recorded in DRG neuron cultures before [12, 19].Besides, they have been observed in many other chemical and neuronal systems [20, 21], and are thereforeof broader interest. We elaborate on the mechanisms of onset and disappearance of MMOs, and comparethem to other extensively analyzed MMO-generating systems.This paper is organized as follows. In Sec. 2, we describe the model and show various patterns of be-haviour using dynamic simulation for selected parameter values. Sec. 3 identifies different parameter regionscorresponding to the different patterns of behavior using one- and two-parameter continuation with externalcurrent as the primary bifurcation parameter and maximal conductance of the Na v In this paper, we focus on a single compartment minimal conductance model. Following [10], the modelaccounts for currents due to two sodium channels: Na v I . ) and Na v I . ), two potassium channels:a delayed rectifier ( I K ) and an A-type transient ( I KA ) channel, and a leak channel ( I l ). These are theprimary ion channels found on the membrane of a small DRG neuron. The equation for membrane voltagedynamics are written in the following Hodgkin-Huxley [22] type of form: C d V d t = I ext A − ( I . + I . + I K + I KA + I l ) , (1)2able 1: Model parameter valuesParameter Value Units A (area) 2168.00 µ m C µ F / cm E Na mV E K -84.70 mV E l -58.91 mV g . mS / cm g . mS / cm g K mS / cm g KA mS / cm g l mS / cm where V is the neuron membrane voltage ( mV ), C is the specific membrane capacitance ( µ F / cm ), t istime ( ms , milliseconds), I ext is the external applied current, and A is the area. I ext /A in Eq. (1) has thedimension µ A / cm . The other specific currents on the right hand side of Eq. (1) are calculated as follows: I . = g . m . h . s . ( V − E Na ) , I . = g . m . h . ( V − E Na ) , I K = g K n K ( V − E K ) , I KA = g KA n KA h KA ( V − E K ) , and I l = g l ( V − E l ) . Therein, g i and E j are the specific maximalconductances ( mS / cm ) and equilibrium potentials ( mV ), respectively, for i = 1 . , . , K,KA, l , and j = N a, K, l .All the activation and inactivation variables x ( x = m . , h . , s . , m . , h . , n K , n KA , h KA ) are di-mensionless variables that can vary between 0 to 1. They are calculated from a corresponding differentialequation of the following form: d x d t = x ∞ ( V ) − xτ x ( V ) , (2)where the nonlinear expressions for x ∞ and τ x are given in the Appendix A.Leak current kinetics, area, membrane capacitance, and equilibrium potential values for a small DRGneuron were extracted from [10]. Maximal conductances g K and g KA were estimated to ensure that theircorresponding currents were 6 nA and 1 nA at 0 mV when the cell is initially depolarized to -120 mV , and g . was set to 18 mS / cm , based on [10]. g . was set to 7 mS / cm , which led to the generation of oneaction potential (current threshold) at 100 pA when the neuron is at the resting membrane potential (RMP),which is determined by simulating the model for I ext = 0 pA . The current threshold of 100 pA was chosenbased on approximate values from previous experiments and simulations [9, 19]. The parameter values ofthis model are listed in Table 1. These parameter values result in an RMP of -66.48 mV which belongs tothe physiological range of RMP recorded in small DRG neurons in [23], and the resulting action potentialamplitude (approximately 120 mV ) is comparable to that reported in [24].The final equation (1) reads: C d V d t = I ext A − ( g . m . h . s . ( V − E Na )+ g . m . h . ( V − E Na )+ g K n K ( V − E K )+ g KA n KA h KA ( V − E K )+ g l ( V − E l )) g . , MMOs are observed. a.:Dynamic simulations for g . at 7 mS / cm , and I ext = 100, 106, 120 pA . b.: Dynamic simulations for g . at 4.5 mS / cm , and I ext = 115, 215, 230 pA . No MMOs are observed in this case.Numerical integration and bifurcation analysis were primarily done in XPPAUT [25] and cross checkedwith MATCONT [26]. In XPPAUT, default settings were used except for the following: NTST = 100,
Method = Stiff,
Tolerance = 1e-7,
EPSL , EPSU , EPSS = 1e-7,
ITMX , ITNW = 20,
PARMIN = 0,
PARMAX = 300. In MATCONT, the following settings were kept:
MaxCorrIters = 20,
MaxTestIters = 20,
FunTolerance = 1e-6,
VarTolerance = 1e-7,
TestTolerance = 1e-7,
NTST = 300, tolerance = 1e-4,
MaxStepsize = 1 for steady state continuation and 10 for periodic solution continuation. Integra-tion was performed using ode15s. Integration option
RelTol was set to 1e-8.In a first step, we present selected dynamic simulations of the above equations to illustrate some char-acteristic patterns of behavior, to be analyzed in more detail in the remainder. Results are shown in Fig. 1.Initial condition for all simulations was the stable steady state for I ext = 0 pA .In the first row of Fig. 1 the maximum conductance of the Na v mS / cm and thevalue of the external current is increased from 100 pA in the left diagram, to 106 pA in the middle, to 120 pA in the right diagram. In the left diagram, for the lowest value of I ext , a stable steady state is attained after thefiring of an action potential, whereas periodic firing of large amplitude action potentials is observed for thehighest value of I ext in the right diagram. For values of I ext in between, there is a region where mixed modeoscillations (MMOs) are observed as illustrated in the middle diagram. There, after some initial transient, aperiodic regime is attained with one large amplitude action potential and eight small amplitude subthresholdpeaks per period. 4 different situation is seen in the second row of Fig. 1 with some representative simulations. There,the maximum conductance of the Na v mS / cm . Again the values of the externalcurrent are increasing from the left to the right. The qualitative behavior in the left and the right diagrams,for the lowest and the highest value of I ext , is similar to the behavior shown in the corresponding diagramsof Fig. 1a. However, in contrast to Fig. 1a, no MMOs are found in the intermediate range of I ext . Weillustrate this in the middle diagram for one specific value of I ext of 215 pA , where the cell potential decaysto a stable steady state after firing of three action potentials. As we will show in more detail in the nextsection, MMOs do not exist for any value of the intermediate range for g . equal to 4.5 mS / cm .Both of the cases shown in Fig. 1 have been observed in DRG culture recordings. See [23] for recordingsresembling Fig. 1b, and [12] and [19] for recordings displaying MMOs as in Fig. 1a. In order to explain the transitions between different dynamical patterns observed in this model, we performone-parameter and two-parameter continuation of solutions, with I ext as the primary bifurcation parameter,and g . as the secondary bifurcation parameter. First, we perform one-parameter continuations of steadystate and periodic solutions upon varying the primary bifurcation parameter I ext . Results are shown in Fig. 2for four different values of g . . The first diagram in Fig. 2a is for a value of 4.5 mS / cm corresponding tothe scenario in Fig. 1b, whereas the third diagram in Fig. 2c is for a value of 7 mS / cm corresponding tothe scenario in Fig. 1a. Two additional scenarios for values of 5 and 8 mS / cm are shown in Figs. 2b and2d. In all the four diagrams of Fig. 2, a branch of stable steady states is obtained for low values of I ext starting from the left boundary of the corresponding diagram. It is indicated by the red solid line andcorrespond to the behavior shown in the left diagrams of Figs. 1a, b. The stable steady states becomeunstable at a subcritical Hopf bifurcation point (HB), from where a branch of unstable periodic solutionsemerges indicated by the blue circles in Fig. 2.Furthermore, in all the four diagrams of Fig. 2, a branch of stable periodic solutions is observed for highvalues of I ext at the right boundary of the corresponding diagrams. It is indicated by the green filled circlesand correspond to periodic firing of action potentials as shown in the right diagrams of Figs. 1a, b. In allthe four cases, these branches of stable periodic solutions lose their stability at a cyclic limit point (CLP ),giving rise to a branch of unstable periodic solutions.Qualitative differences in the four diagrams of Fig. 2 occur with respect to the unstable periodic branchesindicated by the blue circles. In the first two diagrams, the unstable periodic solution branches are connectedand show two more cyclic limit points (CLP and CLP ); after the third cyclic limit point (CLP ), they turninto the stable periodic solution branch with periodic firing of action potentials as described above. Incontrast to this, in the last two diagrams, the unstable periodic solution branches are disconnected and CLP disappears. At the end points of these branches, the period of the unstable periodic solutions is increasingrapidly during continuation with XPPAUT and MATCONT, indicating the presence of a period infinitysolution at the end points.Another difference occurs with respect to the unstable steady state branches indicated by the dashedlines in Fig. 2. They display hysteresis with two limit points (LP and LP ) in Fig. 2c, one of which (LP )has moved out of the positive range for I ext in Fig. 2d.Further, we see that the I ext value of the bifurcation points varies significantly between the four diagramsof Fig. 2, indicating a high sensitivity to g . . For the increasing values of g . from Fig. 2a to 2d, thebifurcation points are shifted to lower values of I ext . Besides the absolute I ext value, we find that the5igure 2: Bifurcation diagrams for g . = 4.5, 5, 7 and 8 mS / cm for diagrams a, b, c and d, respectively.For lower values of g . in diagram (a), MMOs are not observed, and there is a region of bistability betweensteady state and periodic firing of action potentials, as shown by the orange shaded region. This bistability isnot present in diagrams b, c, d. Instead, MMOs are observed in these diagrams in the purple shaded region.MMOs solution branches will be discussed separately in section 3 and are not included in this figure. HB:Hopf bifurcation point, CLP: Cyclic limit point, LP: limit point6elative position of the HB point and the CLP point is of major importance for the qualitative differencesreported in Fig. 1. In Fig. 2a, the I ext value of HB is higher than that of CLP , leading to an overlap betweenstable steady state and stable periodic solutions indicated by the orange region of Fig. 2a. An increase of I ext will lead to the periodic firing of action potentials when the HB point is crossed as shown in the scenario inFig. 1b. A transition back to stable steady states will occur at the value of the CLP point if I ext is decreasedagain afterwards. Between the CLP and the HB point in Fig. 2a, the system is bistable, i.e., the initialconditions regulate whether a stable periodic or a stable steady state solution is attained.The situation is fundamentally different in Figs. 2c and 2d. Here, the I ext value of the HB point is lowerthan the value of the CLP point, leading to a situation where no stable attracting solutions are shown inthe purple shaded region of these diagrams. This is the region where various types of stable periodic andaperiodic MMOs exist. The MMOs solution branches are missing in Fig. 2 and will be discussed in the nextsection. Further, we will show that in these cases, the CLP point provides a strict upper limit of the MMOsregion, whereas the lower limit is not determined by the HB point but by a value close by where the timebetween subsequent large amplitude action potentials of the MMOs tends to infinity.To map out the regions in the I ext and g . parameter space with different patterns of behavior, weperform a two-parameter continuation of the relevant critical points HB, CLP , LP , and LP . The resultsare shown in Fig. 3. As mentioned above, the upper boundary of the MMOs region is the curve of theCLP points, whereas the lower boundary is a solution where the time between subsequent action potentialstends to infinity close to the HB curve. These boundaries are seen best in Fig. 3b. A direct calculation andcontinuation of the lower boundary is substantially challenging and was not done. Instead, we determinethe lower boundary of the MMOs region by point-wise dynamic simulation over a prolonged time period.In summary, we find that the transition from the stable steady state region (‘no pain’) to the repetitive firingof action potentials (‘pain’) differs depending on the value of g . . For high values of g . , MMOs occur.As we will show in the next section, the frequency of action potentials in this region is increasing step bystep as the stimulus I ext is increased. In contrast to this, for values of g . below the intersection of the HBand the CLP curve, we have the orange bistable region with a ‘hard’ onset of the periodic firing of actionpotentials with high frequency.This analysis suggests that the small DRG neuron dynamics depend strongly on the expression ofNa v v g i ( i = 1 . , K, KA ) using two-parameter continuations of the critical bifurcation points with I ext as the primary bifurcation parameterand the corresponding maximal conductance as the secondary bifurcation parameter. Results are shown inFig. 4. As seen in Fig. 4a and Fig. 4c, g . and g KA have negligible effect on the HB and the LP points. TheCLP point is sensitive only to lower values of g KA . All the bifurcation points are sensitive to g K , as seenin Fig. 4b. Here, the CLP , HB, and LP points vary substantially from 0 to 300 pA in a small range of g K . In this section, our focus is on the periodic and aperiodic MMOs solutions already mentioned in the previoussections. For the characterization of periodic MMOs solutions, we apply the nomenclature introduced, forexample, in [27]. Basic MMOs patterns consist of L large amplitude peaks (action potentials) followed by S small amplitude (subthreshold) peaks per period, termed as L S patterns in this notation. In particular, L is equal to one in the remainder of this section. More complex patterns arise due to the concatenationof different basic patterns, for example, a pattern of the form L S L S can occur between the basic patterns7igure 3: Two parameter plot with g . as the secondary continuation parameter. a.: Variation over a largeinterval of g . . b.: Zoomed in version of a. near the intersection of the HB point and the CLP point.Figure 4: Two parameter plot with the following secondary continuation parameters: a.: g . , b.: g K and c.: g KA . 8igure 5: Basic MMOs solutions of the type: a.: , b.: , and c.: for selected values of I ext . Upperrow: temporal evolution of membrane voltage, lower row: orbits in the V , h ka , n K phase space. L S and L S . It consists of L action potentials, followed by S subthreshold peaks, followed by L actionpotentials, followed by S subthreshold peaks in each period.We show some characteristic basic patterns of MMOs in Fig. 5 for different values of I ext . In theremainder of this section, we use the default parameter values from Table 1 and the value of g . is equal to7 mS / cm corresponding to Fig. 2c. According to this figure and our previous results, we expect MMOsin the range of I ext of roughly 103 to 117 pA . More precise values will be given in the course of thisdiscussion. According to the nomenclature mentioned above, the patterns in Fig. 5 can be characterizedas for I ext = 107 pA , for I ext = 110 pA , and for I ext = 114 pA . In this series, the number ofsmall subthreshold peaks is decreasing with increasing external current, and the distance between the actionpotentials is decreasing with increasing external current.For the dynamic simulation of MMOs, it is important to note that the system has multiple time scales.The s . variable is by far the slowest variable. Therefore, we perform all the dynamic simulations in thissection with a startup phase of 100,000 ms to achieve the desired asymptotic behavior of all variables. Thetime window shown, for example, in Fig. 5, starts after this startup phase of 100,000 ms .To add more details to the picture presented in Fig. 5, we perform a one-parameter continuation of thebasic MMOs patterns illustrated in Fig. 5. Results are shown in Fig. 6. The maximum amplitude of theseperiodic solutions is almost constant and therefore not compelling; instead of the amplitude, we use the9igure 6: Basic periodic solution branches with one action potential per period in the range of I ext from 105to 120 pA . Solid lines: stable periodic solutions, dashed line: unstable periodic solutions.period of different solutions for graphical representation of the results.Fig. 6 shows a sequence of isolated periodic solution branches, with one action potential per period. Thenumber of small amplitude peaks between the action potentials and the period increases from the right to theleft in the direction of decreasing external current. The right most branch with label corresponds to theperiodic firing of action potentials without any small amplitude peaks in between, as illustrated in the rightdiagram of Fig. 1a. This periodic solution branch becomes unstable at a cyclic limit point at I ext = 116 . , which corresponds to the CLP point in Fig. 2c. On every other branch in Fig. 6, the correspondingperiodic solution becomes unstable at a cyclic limit point on the left and at a period doubling bifurcationpoint on the right. The values of I ext at the cyclic limit points are indicated by the red lines and the valuesat the period doubling bifurcation points by the green lines in Fig. 6. These values are also listed in Table 2.Solutions below 105 pA are not shown in this figure. Below 105 pA , the distance between the largeamplitude action potentials becomes larger and larger as we increase I ext and finally tends to infinity close10able 2: Values of I ext at the cyclic limit points (CLP) and the period doubling bifurcation points (PD) inFig. 6.Solution type CLP ( pA ) PD ( pA ) I ext equals 102.9935 pA . Accordingly, the number ofsubthreshold peaks becomes larger and larger and their amplitude smaller and smaller as we approach thecritical point. As an example, we show dynamic simulations for a value of I ext = 102 .
992 pA which isslightly below the subcritical Hopf bifurcation point, shown in Fig. 7. The distance between two actionpotentials is roughly 40,000 ms . However, as shown in the phase diagram in Fig. 7b, the orbit is a narrowband and does not seem to be strictly periodic. For I ext slightly below this value, MMOs finally vanish anda stable steady state is obtained.We find concatenated periodic solutions in the gaps of the basic periodic patterns in Fig. 6 between theperiod doubling points and the cyclic limit points of the subsequent solution branches on the right. Wefurther study the dynamic behavior in these regions for selected values of I ext using dynamic simulations.Again, it is crucial to account for the long transient phase introduced by the very slow s . variable asdescribed above. Some characteristic patterns of behavior in the range of 112.9 to 113.2 pA are shown inFig. 8. According to the aforementioned nomenclature, the solution in Fig. 8a can be characterized as aconcatenation between the basic pattern on the left of this value and the basic on the right of this valuein Fig. 6, leading to a solution with two action potentials per period. Accordingly, Fig. 8b demonstratesa (1 ) pattern with 3 action potentials per period, Fig. 8c a (1 ) pattern with 4 action potentials perperiod, and Fig. 8d a (1 ) pattern with 5 action potentials per period.Subsequently, we order the MMOs solutions which were found for selected values of I ext in a tree likestructure in Fig. 9 containing basic and concatenated MMOs patterns as described above. The correspondingvalues of I ext in pA are given in parentheses. The solutions highlighted in yellow correspond to those shownin Fig. 8. It is worth noting that the solution tree is not complete, since only selected values of I ext have beenconsidered. For example, we expect that between the solution (1 ) at I ext = 107.27 pA and the solution (1 ) at I ext = 107.3 pA , another solution of the form (1 ) can be found for some suitable value of I ext , so that the solutions form a regular so-called Farey sequence [27]. Furthermore, we expect that evenhigher order concatenated solutions can be found for some suitable values of I ext .Close to the cyclic limit point at I ext = 116.9811 pA in Fig. 6 corresponding to CLP in Fig. 2c beforethe MMOs disappear, the solution consists of one small amplitude peak and multiple large amplitude peaks.If n is the number of large amplitude peaks, this can be written as a concatenation of one solution and11igure 7: a.: MMOs for I ext = 102.992 pA below the Hopf bifurcation point at I ext = 102.9935 pA . b.:Representation of the solution in the V, h KA , n K phase diagram. ( n −
1) 1 solutions as (1 ) n − . Selected solutions for this region are shown in Fig. 10. The number oflarge amplitude action potentials per period is increasing in this sequence from the left to the right.The solution tree in Fig. 10 is also not complete. For example, we expect that one can also find solutionsof the form (1 ) and (1 ) between (1 ) and (1 ) for some suitable value of I ext in between.For an additional characterization of periodic MMOs, we introduce a firing number F . Following [27], F is defined as the number of small amplitude subthreshold peaks per total number of peaks in a period. Fora basic L S pattern, F is given by: F = SL + S . (3)Accordingly, − F is the firing rate of action potentials per period, which is even more interesting fromthe physiological point of view.The firing number of concatenated MMOs solutions can be calculated using the Farey arithmetic [28].According to this arithmetic, the Farey sum ⊕ of two rational numbers p /q and p /q is defined as: p q ⊕ p q = p + p q + q . (4)Using this definition, the firing number F of a concatenated solution L S L S is, for example, obtainedfrom the following: F = F ⊕ F = S + S L + S + L + S . (5)Furthermore, for the firing numbers of two adjacent solutions in a regular Farey sequence p /q and p /q , the following condition holds: | p q − p q | = 1 . (6)An illustration of the Farey arithmetic for specific MMOs solutions sequence is shown in Table 3.12igure 8: A sequence of concatenated periodic solutions. a.: at I ext = 112 . , b.: (1 ) at I ext = 113 . , c.: (1 ) at I ext = 113 .
18 pA , d.: (1 ) at I ext = 113 . .13 (105.64) (106) (105.85) ( ) (105.87) ( ) (105.88) (106.6) (106.2) ( ) (106.23) (106.8) (106.65) ( ) (106.69) ( ) (106.7) (106.8) (107.4) (107.2) ( ) (107.27) ( ) (107.3) (108.2) (108) ( ) (108.05) (109.3) (109.1) ( ) (109.15) (109.3) (110.9) (110.6) ( ) (110.7) (113.3) (112.9) ( ) (113.1) ( ) (113.18) ( ) (113.2) (116.99) (116.6) ( ) (116.8) ( ) (116.85)( ) (116) Figure 9: Tree of selected periodic MMOs solutions. Numbers in parentheses are values of I ext in pA corresponding to the solution on top of it. Solutions highlighted in yellow are shown in Fig. 8.14 (113.3) (116.99) (116.6) ( ) (116.8) ( ) (116.85) ( ) (116.9) ( ) (116.95) ( ) (116.965) ( ) (116.97) ( ) (116.972) ( ) (116.975) Figure 10: Selected periodic MMOs patterns observed below but close to the cyclic limit point CLP inFig. 2c before small amplitude oscillations disappear. Numbers in parentheses are the corresponding valuesof I ext in pA , corresponding to the solution on top of it.Table 3: An illustration of MMOs solution sequences satisfying the Farey arithmetic.MMOs solution Firing number p q − p q ⊕ (1 ) ⊕ (1 ) ⊕ I ext slightly above the period doubling points. This is illustrated in Fig. 11 by two simulations.The diagrams on the left demonstrate the dynamic behavior after the startup phase at I ext = 108.9 pA , belowthe period doubling point at I ext = 108.9962 pA , with a stable periodic MMOs solution. The diagramson the right demonstrate a second solution at I ext = 109 pA after the startup phase, slightly above the perioddoubling point. The aperiodic, seemingly chaotic behavior, is not so obvious from the voltage dynamics;however, irregularity is seen in the dynamics of the s . variable. In Fig. 11c, s . forms a thick straight bandover a long time period, implying that no variation is seen in the s . oscillations. However, in Fig. 11d,irregularity in s . is found over this long time period and no repeating patterns are observed. In this work, we attempt to understand the dynamics of a 9D model representative of a small DRG neuron.Small DRG neurons are primary nociceptors and can sense pain. Any damage to them due to injuries, dis-eases, or genetic disorders, can lead to conditions such as loss or gain of nociceptive pain sensation, andneuropathic or inflammatory pain. A bifurcation analysis of this model can aid in understanding the transi-tion of this system from steady state to mixed-mode oscillations, and finally to full blown periodic firing ofaction potentials, where oscillations of any frequency indicate pain of a specific form and intensity [7, 8].The model displays rich dynamics, which we investigate by studying the bifurcations numerically us-ing the external applied current as the primary bifurcation parameter and the maximal conductances g i ( i = 1 . , . , K, KA ) of the sodium and potassium channels as secondary parameters. We show that, inparticular, g . and g K are the most sensitive maximal conductances. We provide a detailed analysis for g . as the secondary parameter. We show that there is a hard onset of periodic firing of action potentials dueto hysteresis between stable steady state and periodic firing of action potentials for low values of g . . Thispattern of behavior can also be found in the original Hodgkin-Huxley equations (see, for example, [3]). Forhigh values of g . , the frequency of firing of action potentials is increasing step by step as we pass througha region of MMOs where the distance between the action potentials is getting smaller and smaller as thenumber of subthreshold peaks between the action potentials is reduced step by step until they finally vanish.Although the region of MMOs is rather narrow in terms of I ext for the parameter values considered in thispaper, it represents a second, and a fundamentally different path to pain.Using selected dynamic simulations, we conjecture that the periodic MMOs build Farey sequences. SuchFarey sequences have also been observed for various other systems displaying MMOs (see [27, 29–32] forexamples); however, they have not been widely studied for neuron models (see [33] for an example). Giventhe diversity and abundance of neuron models that can generate MMOs [20], it will be interesting to explorethe existence of Farey sequences in other neuron models as well.Besides periodic, we also find aperiodic MMOs for very small ranges of I ext . We conjecture that theseaperiodic MMOs solutions are chaotic. Further investigations are required to validate this hypothesis. Froma mechanistic point of view, it would be interesting to record such chaotic behavior in DRG cultures andfind its implications on pain sensation.From the mathematical point of view, the 9D model used in this study is rather complex and prohibitsfurther analytical insight as demonstrated for example in [34] for a lower dimensional problem. To gainfurther theoretical insight it would therefore be desirable to reduce the present model to a lower dimensionalproblem showing similar patterns of behavior.From the physiological point of view, the model used in this study is still relatively simple. Towards amore realistic description of small DRG neurons, additional ion channels should be taken into account such16igure 11: Simulations before and after the period doubling bifurcation at I ext = 108.9962 pA . Left column: I ext = 108 . , right column: I ext = 109 pA . After the period doubling bifurcation, the system exhibitschaotic-like behavior which is evident in the dynamics of s . .17s Na v I ext = 0 pA ), others do not [24]. This indicates thatsimply fixing maximal conductances as constants will not suffice, and there is a need to analyze an ensembleof possible parameter values to capture the heterogeneity observed in electrical recordings.Mathematical understanding of the sensing of pain is necessarily an evolving process of manipulatingmodel scale to suitably match the minimum physiological details associated with pain. Each step in thisprocess involves comparing predictions with experimental observations and identifying how parametersconnected with various ion channels relate to pain. Thus both model elaboration and reduction are potentialfuture areas of interest, and can enable a rigorous investigation of possible dynamics that can be displayedby this system. This can further shape our understanding of pain sensation and how it can be controlled. A Appendix
A.1 Na v α m . = 15 .
51 + exp (cid:18) ( V − − . (cid:19) (7) β m . = 35 .
21 + exp (cid:18) V + 72 . . (cid:19) (8) α h . = 0 . (cid:18) V + 122 . . (cid:19) (9) β h . = − . . (cid:18) ( V + 5 . − . (cid:19) (10) α s . = 0 . . (cid:18) V + 93 . . (cid:19) (11) β s . = 132 . − .
051 + exp (cid:18) V − . . (cid:19) (12)For x = m . , h . , s . : x ∞ ( V ) = α x ( V ) α x ( V ) + β x ( V ) , (13)and τ x ( V ) = 1 α x ( V ) + β x ( V ) (14)18he kinetics of Na v m . corresponds to the activation gating variable, h . to the fast-inactivation gating variable and s . to the slow-activation gating variable. A.2 Na v α m . = 2 . − . (cid:18) V − . . (cid:19) (15) β m . = 7 . (cid:16) V + 46 . . (cid:17) (16)For x = m . : x ∞ ( V ) = α x ( V ) α x ( V ) + β x ( V ) , (17)and τ x ( V ) = 1 α x ( V ) + β x ( V ) (18) τ h . = 1 .
218 + 42 . × exp (cid:18) − ( V + 38 . × . (cid:19) (19) h . ∞ = 11 + exp (cid:18) V + 32 . (cid:19) (20)The kinetics of Na v m . and h . are similar activation and inactivationgating variables, respectively. A.3 K equations α n K = 0 . × ( V + 14 . − exp (cid:18) V + 14 . − (cid:19) (21)with α n K = 0 . × for V = − . . β n K = 0 . × exp (cid:18) V + 55 − . (cid:19) (22) n K ∞ = 11 + exp (cid:18) − ( V + 14 . . (cid:19) (23) τ n K = 1 α n K + β n K + 1 (24)The kinetics of K channel were taken from [35]. 19 .4 KA equations n KA ∞ =
11 + exp (cid:18) − ( V + 5 . . (cid:19) (25) τ n KA = 0 .
25 + 10 . × exp (cid:18) − ( V + 24 . × . (cid:19) (26) h KA ∞ = 11 + exp (cid:18) V + 49 . . (cid:19) (27) τ h KA = 20 + 50 × exp (cid:18) − ( V + 40) × (cid:19) (28)The kinetics of KA channel were taken from [9] Acknowledgements
This project was supported, in part, with support from the Indiana Clinical and Translational Sciences In-stitute funded, in part by Award Number UL1TR002529 from the National Institutes of Health, NationalCenter for Advancing Translational Sciences, Clinical and Translational Sciences Award. The content issolely the responsibility of the authors and does not necessarily represent the official views of the NationalInstitutes of Health. The authors also thank Dr. Haroon Anwar, New Jersey Institute of Technology, USA,for helping with model selection and building, and for reviewing this manuscript; Max Planck Institute forDynamics of Complex Technical Systems, Magdeburg, Germany, for sponsoring trips to strengthen the col-laboration; and Muriel Eaton and Dr. Yang Yang, Purdue University, USA, for insightful discussions onDRG neurons and pain sensation.
Conflict of interest
The authors declare that they have no conflict of interest.
References [1] William C Troy. The bifurcation of periodic solutions in the Hodgkin-Huxley equations.
Quarterly ofApplied Mathematics , 36(1):73–83, 1978.[2] Brian Hassard. Bifurcation of periodic solutions of the hodgkin-huxley model for the squid giant axon.
Journal of Theoretical Biology , 71(3):401–420, 1978.[3] John Rinzel. Numerical calculation of stable and unstable periodic solutions to the Hodgkin-Huxleyequations.
Mathematical Biosciences , 49:27–59, 1980.[4] Eugene M Izhikevich.
Dynamical systems in neuroscience . MIT press, 2007.205] Huiwen Ju, Alexander B. Neiman, and Andrey L. Shilnikov. Bottom-up approach to torus bifurcationin neuron models.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 28(10):106317, 2018.[6] Charles S Sherrington. Qualitative difference of spinal reflex corresponding with qualitative differenceof cutaneous stimulus.
The Journal of physiology , 30(1):39, 1903.[7] Laiche Djouhri, Stella Koutsikou, Xin Fang, Simon McMullan, and Sally N. Lawson. Spontaneouspain, both neuropathic and inflammatory, is related to frequency of spontaneous firing in intact c-fibernociceptors.
Journal of Neuroscience , 26(4):1281–1292, 2006.[8] Adrienne E. Dubin and Ardem Patapoutian. Nociceptors: the sensors of the pain pathway.
The Journalof Clinical Investigation , 120(11):3760–3772, 11 2010.[9] Patrick L Sheets, James O Jackson, Stephen G Waxman, Sulayman D Dib-Hajj, and Theodore R Cum-mins. A Nav1. 7 channel mutation associated with hereditary erythromelalgia contributes to neuronalhyperexcitability and displays reduced lidocaine sensitivity.
The Journal of physiology , 581(3):1019–1031, 2007.[10] Jin-Sung Choi and Stephen G Waxman. Physiological interactions between Nav1. 7 and Nav1. 8sodium channels: a computer simulation study.
Journal of neurophysiology , 106(6):3173–3184, 2011.[11] Danielle Sundt, Nikita Gamper, and David B. Jaffe. Spike propagation through the dorsal root gangliain an unmyelinated sensory neuron: a modeling study.
Journal of Neurophysiology , 114(6):3140–3153,2015. PMID: 26334005.[12] Ron Amir, Martin Michaelis, and Marshall Devor. Membrane potential oscillations in dorsal rootganglion neurons: Role in normal electrogenesis and neuropathic pain.
Journal of Neuroscience ,19(19):8589–8596, 1999.[13] Anthony M. Rush, Theodore R. Cummins, and Stephen G. Waxman. Multiple sodium channels andtheir roles in electrogenesis within dorsal root ganglion neurons.
The Journal of Physiology , 579(1):1–14, 2007.[14] Elliot S. Krames. The Role of the Dorsal Root Ganglion in the Development of Neuropathic Pain.
Pain Medicine , 15(10):1669–1685, 10 2014.[15] Temugin Berta, Yawar Qadri, Ping-Heng Tan, and Ru-Rong Ji. Targeting dorsal root ganglia andprimary sensory neurons for the treatment of chronic pain.
Expert opinion on therapeutic targets ,21(7):695–703, 2017.[16] Young-Ah Rho and Steven A. Prescott. Identification of molecular pathologies sufficient to causeneuropathic excitability in primary somatosensory afferents using dynamical systems theory.
PLOSComputational Biology , 8(5):1–14, 05 2012.[17] St´ephanie Ratt´e, Yi Zhu, Kwan Yeop Lee, and Steven A Prescott. Criticality and degeneracy in injury-induced changes in primary afferent excitability and the implications for neuropathic pain.
Elife , 3,2014.[18] Parul Verma, Achim Kienle, Dietrich Flockerzi, and Doraiswami Ramkrishna. Using bifurcation the-ory for exploring pain.
Industrial & Engineering Chemistry Research , 0(0):null, 0.2119] Qin Zheng, Dong Fang, Jie Cai, You Wan, Ji-Sheng Han, and Guo-Gang Xing. Enhanced excitabilityof small dorsal root ganglion neurons in rats with bone cancer pain.
Molecular pain , 8(1):24, 2012.[20] Morten Brøns, Tasso J. Kaper, and Horacio G. Rotstein. Introduction to focus issue: Mixed modeoscillations: Experiment, computation, and analysis.
Chaos: An Interdisciplinary Journal of NonlinearScience , 18(1):015101, 2008.[21] Mathieu Desroches, John Guckenheimer, Bernd Krauskopf, Christian Kuehn, Hinke M Osinga, andMartin Wechselberger. Mixed-mode oscillations with multiple time scales.
Siam Review , 54(2):211–288, 2012.[22] A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its applicationto conduction and excitation in nerve.
The Journal of Physiology , 117(4):500–544, 1952.[23] Jianying Huang, Malgorzata A Mis, Brian Tanaka, Talia Adi, Mark Estacion, Shujun Liu, SuellenWalker, Sulayman D Dib-Hajj, and Stephen G Waxman. Atypical changes in DRG neuron excitabil-ity and complex pain phenotype associated with a Na v 1.7 mutation that massively hyperpolarizesactivation.
Scientific reports , 8(1):1811, 2018.[24] Yang Yang, Jianying Huang, Malgorzata A. Mis, Mark Estacion, Lawrence Macala, Palak Shah,Betsy R. Schulman, Daniel B. Horton, Sulayman D. Dib-Hajj, and Stephen G. Waxman. Nav1.7-A1632G mutation from a family with inherited erythromelalgia: Enhanced firing of dorsal root ganglianeurons evoked by thermal stimuli.
Journal of Neuroscience , 36(28):7511–7522, 2016.[25] Bard Ermentrout.
Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT forresearchers and students , volume 14. Siam, 2002.[26] A. Dhooge, W. Govaerts, Yu. A. Kuznetsov, H. G.E. Meijer, and B. Sautois. New features of the soft-ware MatCont for bifurcation analysis of dynamical systems.
Mathematical and Computer Modellingof Dynamical Systems , 14(2):147–175, 2008.[27] J. Masełko and Harry L. Swinney. Complex periodic oscillations and Farey arithmetic in the Be-lousovZhabotinskii reaction.
The Journal of Chemical Physics , 85(11):6430–6441, 1986.[28] Godfrey Harold Hardy, Edward Maitland Wright, et al.
An introduction to the theory of numbers .Oxford university press, 1979.[29] F. N. Albahadily, John Ringland, and Mark Schell. Mixedmode oscillations in an electrochemicalsystem. I. A Farey sequence which does not occur on a torus.
The Journal of Chemical Physics ,90(2):813–821, 1989.[30] T Hauck and FW Schneider. Chaos in a Farey sequence through period doubling in the peroxidase-oxidase reaction.
The Journal of Physical Chemistry , 98(8):2072–2077, 1994.[31] Andrzej L. Kawczyski and Peter E. Strizhak. Period adding and broken Farey tree sequence of bifur-cations for mixed-mode oscillations and chaos in the simplest three-variable nonlinear system.
TheJournal of Chemical Physics , 112(14):6122–6130, 2000.[32] Wieslaw Marszalek. Circuits with oscillatory hierarchical Farey sequences and fractal properties.
Circuits, Systems, and Signal Processing , 31(4):1279–1296, Aug 2012.2233] Martin Krupa, Nikola Popovi, Nancy Kopell, and Horacio G. Rotstein. Mixed-mode oscillations in athree time-scale model for the dopaminergic neuron.
Chaos: An Interdisciplinary Journal of NonlinearScience , 18(1):015106, 2008.[34] Jonathan Rubin and Martin Wechselberger. Giant squid-hidden canard: The 3d geometry of theHodgkin–Huxley model.
Biological Cybernetics , 97(1):5–32, Jul 2007.[35] JH Schild, JW Clark, M Hay, D Mendelowitz, MC Andresen, and DL Kunze. A-and C-type rat nodosesensory neurons: model interpretations of dynamic discharge characteristics.