Bifurcation analysis on the reduced dopamine neuronal model

1.   Introduction

The activity of midbrain dopamine neurons, as reflected in levels of extracellular dopamine concentration and the BLOD-fMRI (blood oxygen level dependent functional magnetic resonance imaging, BOLD-fMRI) signals in their target areas, is hypothesized to represent a reward prediction error or confidence in a prediction of a desired outcome^[1,2,3]. Firing patterns of dopamine neurons influence dopamine signaling, for example, electrical stimulation of dopaminergic cells with 40 Hz (Hertz) increased extracellular concentrations of dopamine in the rat striatum more effectively than the same amount of stimulation with 10 Hz^[4,5]. In vitro experiments, dopamine neurons exhibit regular firing with frequency from 1 to 7 Hz^[6,7,8,9]. However, in living organisms, different distribution patterns are observed, and these include regular single spiking firing, irregular single spiking firing, and bursting, which are present in both freely moving^[10,11] and anesthetized^[12,13] rats.

Under different conditions, dopamine neurons exhibit a variety of oscillatory patterns. For example, after blocking the spiking with TTX (Tetrodotoxin), the neurons exhibit a slow wave firing that is approximately sinusoidal and calcium-mediated ^[8,14,15]. This oscillation is in the frequency range of the spike. On the other hand, in both in vivo and ex vivo experiments, blocking^[16,17,18] or minimizing^[8,19,20] the small conductance SK (small conductance calcium activated potassium, SK) potassium channels increases the tendency of the bursting to rhythmic firing rhythms. These bursts usually end in a depolarizing blockade, similar to the bursting seen in rats on long-term antipsychotic drugs^[21,22]. Depolarizing blockage refers to a state in which the resting state of the membrane potential is programmed to support a more depolarization of the active potential generation. Clues to the bursting mechanism include: (1) blocking the spike firing and SK channels and triggering a plateau potential oscillation in the bottom layer^[8,23,24], during which the depolarization phase can last for several seconds, and (2) L-type calcium channel agonists are also sufficient to trigger bursting behavior. L-type calcium channel-sparing agents disrupt the bursting behavior triggered by the SK channel blockade^[25,26]. The dopamine neuron model proposed by Yu and Canavier^[27,28] is the only model that captures the behavior of the SK channel blocking to mediate bursting in the dopamine neuron. The most notable feature of the model is that the resting phase is more depolarized than the membrane potential during the inter-spike.

In many mathematical models of bursting emanation, the membrane potential has a slow oscillation below, during which the spike of the bursts appear during depolarization, and the inter-spike intervals without emanation are more hyperpolarized than the inter-spike interval, which are often called square-wave bursts^{[29,30,31,32]}. A class of anomalous bursting is considered, in which the non-emitting inter-spike intervals are more depolarized than the mean membrane potential observed during the interspike interval. The single atrial compartment model of dopamine neurons studied in this chapter is derived from the model constructed by Yu and Canavier^[33,34].

Due to the complexity of neurons, neuron models are often characterized by high-dimensional complex nonlinear dynamical systems. The research process of a high-dimensional model is very difficult, so we use the method in the literature^[35]: a projection method to reduce dimension. The dopamine neuron model studied in this chapter is a 13-dimensional nonlinear dynamical system. The original 3D system simulates the biological properties of the dopamine neuron very well. However, it is very difficult to study the kinetic properties of this model directly. We first use variable correlation to reduce the original 3D model in two steps to obtain a simplified system. Then, we discuss the variation rule of the spikes in the simplified system. By changing the parameters $V_{Leak}$ and $g_{Na}$ simultaneously, it is found that the number of spikes in the simplified model of the action potential varies regularly during a single week, and a number of turbid regions are created. Finally, the birfucation of equilibrium points of the simplified system are studied. The remaining two-dimensional bifurcation of the simplified system, such as the CP bifurcation (cusp bifurcation), the BT bifurcation (Bogdanov-Takens bifurcation), and the GH bifurcation (generalized Hopf bifurcation), are obtained by calculation. A detailed analysis of the topological structure in the vicinity of the BT point is carried out to obtain the specific structure of the bifurcation on the central manifold in the vicinity of the point.

2.   Materials and methods

The model presented in this paper is a single-compartment midbrain dopamine neuron (cylinder of diameter $d = 15\; \mathtt{µ}\mathrm{m}$ and length $L = 25\; \mathtt{µ}\mathrm{m}$) described by the following transmembrane potential^[33]:

where $I_{Na} = g_{Na}m^3hh_s(V-V_{Na})$ and $I_{CaL} = g_{CaL}l(V-V_{Ca})$ are fast sodium and L-type calcium currents, and $I_{KDR} = g_{KDR}n^3(V-V_K)$, $I_{KA} = g_{KA}p(q_1/2+q_2/2)(V-V_K)$, $I_{KERG} = g_{KERG}o(V-V_K)$, $I_{KSK} = g_{KSK}(V-V_K)/(1+(0.00019/[Ca])^4)$ are delayed rectifier, A-type, ether-a-go-go-related, and small conductance calcium dependent potassium currents, respectively. $I_H = g_Hm_H^2(V-V_H)$ is a nonspecific hyperpolarization-activated cation current and leak current $I_{Leak}$ is comprised of a calcium ion specific component $I_{LCa} = g_{LCa}(V-V_{Ca})$ and a nonspecific component $I_{LNS} = g_{LNS}(V-V_{Leak})$. $I_s$ denotes a small external stimulus current. In this paper, we assume that without the external stimuli, namely, $I_s = 0$, the membrane capacitance $C_m = 1\; \mathtt{µ} \mathrm{F}/\mathrm{cm^2}$.

The voltage-dependent conductances are described using the Hodgkin-Huxley formalation:

with $x_\infty = 1/[1+\exp(-\frac{V-x_{half}}{x_k})]$, in which $x_{half}$ and $x_k$ represent the half-activation voltage and slope of the steady-state Boltzmann fit to $x_\infty$. Time constants:

The current $I_{KERG}$ uses a kinetic scheme described previously^[33]. The current description requires following two differential equations:

where $o$ and $i$ represent the fraction of ERG (Ether-a-go-go-RelatedGene) channels in the open and inactivated states, and $\alpha_o$, $\beta_o$, $\alpha_i$, $\beta_i$ denote the corresponding reaction rates:

The $Ca^{2+}$ balance is described by:

where $[Ca]$ denotes the $Ca^{2+}$ concentration in mM (mM/L), and $f_{Ca} = 0.018$ and $F$ are the fraction of the unbuffered free calcium and Faraday constant. $I_{Cap}$ represents extrusion of $Ca^{2+}$, which is modeled by a non-electrogenic pump:

where $I_{Capmax} = 11\; \mathtt{µ} \mathrm{A}/\mathrm{cm^2}$.

The parameters for each current are listed in Table 1.

Table 1.  Parameter values for the dopamine neuron model.

Maximal conductances ($\mu S/cm^2$)	$g_{CaL}=139$, $g_{KDR}=1117$, $g_{KA}=1680$, $g_{KSK}=70$, $g_{Na}=6000$, $g_{KERG}=130$, $g_H = 78$, $g_{LNS} = 280$, $g_{LCa} = 2.45$
Reversal potentials (mV)	$V_{Na}=60$, $V_{Ca}=50$, $V_K=-90$, $V_H=-29$, $V_{Leak}=-65$
Half-activation voltages (mV)	$m_{half}=-30.09$, $h_{half}=-54$, ${h_s}_{half}=-54.8$, $n_{half}=-25$, $l_{half}=-45$, ${m_H}_{half}=-77.6$, $p_{half}=-35.1$, ${q_1}_{half}={q_2}_{half}=-80$
Slopes (mV)	$m_k=13.2$, $h_k=-12.8$, ${h_s}_k=-1.57$, $n_k=12$, $l_k=7.5$, ${m_H}_k=-17.317$, $p_k=13.4$, ${q_1}_k={q_2}_k=-6$

---

 | Show Table

DownLoad: CSV

Simulations and dimensionality reduction for the original 13-dimensional single-compartment model was performed using the NEURON software^[36], and data processing was conducted using ORIGIN^[37] and MATLAB softwares^[38]. The bifurcation diagrams for the slow plateau potential was calculated with XPPAUT. The bifurcation diagrams for the reduced three-dimensional system were generated by the MATCONT package^[39].

3.   Results

3.1.   Model reduction

Combined Eqs (2.1)–(2.5) obtaining the original 13-dimensional model, using the formulas and parameters in Section 2, simulated the discharge properties of the original model; for specifics, see references^[33].

Given the dimension reduction experience of the Hodgkin-Huxley model, the four-dimensional reduced to three-dimensional, see references^[35,40]. To start, consider the alternative method with the appropriate variable reduced dimensionality.

In order to better hold and explore the firing properties of the original model, this paper adopts two-step dimensionality reduction. The first step is reduced to eight-dimension and keeps the firing characteristics of the original model; in order to be convenient for dynamic analysis of the reduced model, the second step is reduced to three-dimensional.

Step one: Use membrane potential dimensionality reduction. Based on the initial parameters of the original model, draw a phase portrait of each variable and membrane potential (Figure 1). As can be seen from Figure 1, gating variables $m$, $q_1$, $q_2$, $o$, $i$ are closely related to membrane potential, but others have no significant relationship. Utilizing software ORIGIN curve fitting can be drawn:

Figure 1.  Phase portraits of gating variables and membrane potential. Red lines are the corresponding fitted curves ($m$, $q_1$, $q_2$, $o$, $i$).

DownLoad: Full-Size Img PowerPoint

Take these five variables into the original equation reduced to the 8-dimensional model. In order to make discharge properties of the 8-dimensional reduced model consistent with the original model, we adjust $g_{Na} = 5400\; \mu S/cm^2$, $V_K = -85\; mV$ (Millivolt), and other parameters remain unchanged. As can be seen from Figure 2, when spiking, firing properties of the 8-dimensional model and the original model are almost exactly the same (frequency and amplitude). $g_{KSK} = 0$, the original model appears bursting; $g_{KSK} = 17.905\; \mu S/cm^2$, the 8-dimensional model appears bursting, which is consistent with the original model (frequency, period, duty cycle).

Figure 2.  Firing patterns of the original model and 8-dimensional model. (a) spiking of the original model with $g_{KSK} = 0\; \mu S/cm^2$, the frequency is 3 Hz; (b) bursting of the original model with $g_{KSK} = 0\; \mu S/cm^2$, the period is 6.98 $S$, the Duty ratio is 0.65; (c) spiking of the the 8-dimensional model with $g_{KSK} = 17.905 \; \mu S/cm^2$, the frequency is 4 Hz; (d) bursting of the the 8-dimensional model with $g_{KSK} = 17.905 \; \mu S/cm^2$, the period is 6.72 $S$, the Duty ratio is 0.61.

DownLoad: Full-Size Img PowerPoint

Step two: Use gating variables $n$ and $[Ca]$ dimensionality reduction. Similar to the first step, draw a phase portrait of remaining variables and $n$ or $[Ca]$. It can be seen from phase portraits, gating variables $h$, $l$, $p$ strongly depends on $n$, and $h_s$, $m_H$ depends on $[Ca]$. Similarly, it can be obtained by curve fitting:

Take these five variables into the 8-dimensional model and obtain a reduced three-dimensional system as follows:

where specific ion channel models are shown above.

In order to better explore dynamic properties of the system Eq (3.1), adjust parameters $I_{Capmax} = 21\; \mu A/cm^2$, $g_{Na} = 6000\; \mu S/cm^2$ (the original value), $g_{CaL} = 60\; \mu S/cm^2$, $g_{KDR} = 17\; \mu S/cm^2$, and other parameters remain unchanged.

3.2.   Influence of two-parameter on firing pattern

First, we consider the impact of parameter changes on the firing pattern of the three-dimensional system. By examining the system's firing patterns under various parameters, it is found that they are essentially similar. In this paper, we select the parameters $V_{Leak}$ and $g_{Na}$, and observe the changes in the firing pattern of the simplified model as the parameters change. By statistically analyzing the number of peaks within a single cluster for each combination of $V_{Leak}$ and $g_{Na}$ parameters, we obtain the distribution of the number of spikes within a single burst of the system's firing on a two-dimensional plane ($V_ {Leak}, g_ {Na}$) (Figure 3).

Figure 3.  Spike-counting diagram with two-parameter $V_{Leak}$ and $g_{Na}$. The color bar in the right column shows the number of spikes within a burst. The number of 0 signifies no firing; the number 1 represents spiking; the number 2–34 indicates regular bursting, while the number 35 denotes chaotic bursting.

DownLoad: Full-Size Img PowerPoint

Figure 3 shows the transient pattern of the discharge mode of the system with the change of $g_{Na}$ and $V_{Leak}$. In Figure 3, the horizontal and vertical axes are $g_{Na}$ and $V_{Leak}$, with values ranging from −70 to −50 $\mathrm{mV}$ and 0–6000 $\mathtt{µ} \mathrm{S}/\mathrm{cm}^{2}$, respectively. From Figure 3, it can be seen that with the increase of the parameters, the system exhibits a firing pattern from resting to spiking, then to bursting, and finally to resting again. The different colors in Figure 3 represent the number of spikes in a single burst, and the specific number of spikes indicated by each color is shown on the right side of the ribbon marking the value, in which the light blue color marked with 0 represents that the system is not firing, 1 represents that the peaks are spiking, 2–34 represents that the system is in the regular bursting state at this time, and 35 indicates that the system enters into a chaotic bursting state. When $V_{Leak} = -58\; \mathrm{mV}$, with the increase of $g_{Na}$, when $g_{Na} > 1850\; \mathtt{µ} \mathrm{S}/\mathrm{cm}^{2}$, the system shows a regular additive period firing phenomenon, the number of spikes in a single burst starts to increase gradually from 2, and the different numbers of regions are separated by the chaotic region, which means that each time the number of peaks changes, the system has to go through a chaotic region. Therefore, the 'comb' shaped chaotic region can be clearly seen in Figure 3. In addition, the system also generates three chaotic regions around $V_{Leak}$ values of −55 to −54, −52 to −51$\; \mathrm{mV}$, and values of $g_{Na}$ are from 1000 to 2000$\; \mathtt{µ} \mathrm{S}/\mathrm{cm}^{2}$, and there is also a conical chaotic region between values of $V_{Leak}$ from −63 to −55$\; \mathrm{mV}$ and $g_{Na}$ values of 2000 to 6000$\; \mathtt{µ} \mathrm{S}/\mathrm{cm}^{2}$.

Second, we branch out to system Eq (3.1). On the basis of dimensionality reduction, we select the two parameters $V_{Leak}$ and $g_{Na}$ above for the co-dimension-2 bifurcation analysis, and the results are shown in Figure 4. Figure 4(a) shows the co-dimension-2 bifurcation diagrams for the values of $V_{Leak}$ in the range of [−70, 50] $\mathrm{mV}$ and g in the range of [0, 6000] $\mathtt{µ} \mathrm{S}/\mathrm{cm}^{2}$. Figure 4(b) is a localized enlargement of Figure 4(a). From Figure 4(a), it can be seen that the system exhibits a rich variation of equilibrium points as the parameters are varied. From the figure, it can be seen that the co-dimension-2 bifurcation points of the system appear in CP (−57.549678,961.58291), BT (−63.413152, 2476.1036), GH1 (−60.217081, 5015.182), and GH2 (−55.831514,778.22531), respectively.

Figure 4.  Two-parameter $(V_{Leak}, g_{Na})$ bifurcation diagram of the oscillatory regime. (a) Global view; (b) a partial amplification of panel (a). The blue curve, SN, corresponds to the saddle-node bifurcation. The red curve H marks the Hopf bifurcation curve. The green curve PD denotes the different period-doubling bifurcation curves of different limit cycle, the yellow curve HC denotes the homoclinic bifurcation curve, the light blue LPC curve denotes the folding bifurcation curve of the limit cycle, the CP point denotes the Cusp bifurcation, the BT point denotes the Bogdanov-Takens bifurcation, the GH1 and GH2 points denote the Generalized Hopf bifurcation, and the CPC denotes the Cusp bifurcation of the limit cycle.

DownLoad: Full-Size Img PowerPoint

Then, we superimpose Figure 3 on Figure 4 to investigate the effect of system bifurcation on the firing behavior, and the results are shown in Figure 5. It can be seen that the chaotic region $I_{1}$ is formed by the PD bifurcation, because when the system generates PD bifurcation, the PD cascade is formed by it. With the increase of $V_{Leak}$, another chaotic region $I_{2}$ is formed when the system passes through the LPC, and with the decrease of parameter $V_{Leak}$, $I_{3}$ region is formed when the system passes through the H-curve. As can be seen from Figure 5, when the system passes through the Hopf bifurcation curve from resting to spiking, the system changes from a resting state to a firing state, and the system starts the firing behavior; when the system passes through the PD curve, the spiking changes to a bursting state; and when the system passes through the H curve from left to right, the system transitions to either a resting or a spiking state.

Figure 5.  Influence of two-parameter on firing pattern.

DownLoad: Full-Size Img PowerPoint

Figure 6 uses the bifurcation diagram of the sequence of the inter-spike interval to further delineate the change of the firing behavior of the system. Figure 6(a) presents the firing pattern of the system exhibiting a PD bifurcation when $V_{Leak} = -58\; mV$, Figure 6(b) shows the adding-period bifurcation related to $g_{Na}$, and Figure 6(c) shows the maximum Lyapunov exponent. It is seen that the system exhibits the formation of chaotic region from the PD bifurcation with the increase of g, which corresponds well to the maximum Lyapunov exponent. Figure 6(d) is a local enlargement of Figure 6(a), which shows a clearer bifurcation of the PD bifurcation phenomenon.

Figure 6.  ISI (interspike interval) bifurcation diagram versus $g_{Na}$. (a) The PD bifurcation diagram; (b) adding PD diagram; (c) the largest Lyapunov exponent corresponding PD bifurcation diagram in panel (a); (d) the partial enlargement of the PD bifurcation diagram in panel(a).

DownLoad: Full-Size Img PowerPoint

3.3.   Bifurcation analysis

The Bogdanov-Takens bifurcation point $BT$ appears when $(V_{Leak}, g_{Na})^T = (-63.413152, 2476.1036)^T\triangleq\mu_0$, with corresponding equilibrium coordinate $(V, n, [Ca])^T = (-56.177560, 0.069258865, 0.000022315168)^T\triangleq X_0$. Rewrite the system Eq (3.1) as:

where $X = (V, n, [Ca])^T$, $\mu = (V_{Leak}, g_{Na})^T$, and

Consider the Taylor series of $F(X, \mu)$ around $(X_0, \mu_0)$:

Note matrix

The matrix $A$ has three eigenvalues, that is, $0$, $0$, and $-0.0363794162$. Let $P = (p_1, p_2, p_0)$ be an invertible matrix, which satisfies $P^{-1}AP = J$, where $J = \left(\begin{array}{cc} J_0 & 0 \\ 0 & J_1 \\ \end{array} \right)$, $J_1 = -0.0363794162$, $p_1$ and $p_2$ are generalized eigenvectors of the matrix $A$ corresponding to the double-zero eigenvalues, and $p_0$ contains the generalized eigenvectors of the matrix $J_1$. We have

If $P^{-1} = (q_1, q_2, q_0^T)^T$, then

By calculating expressions (28) and (29) in ^[41], then we obtain

If we choose $\lambda_1$ and $\lambda_2$ as bifurcation parameters, where $\lambda_1 = V_{Leak} + 63.413152$, $\lambda_2 = g_{Na} - 2476.1036$, then

Using the Theorem 1 in ^[41], the reduced system (6) near ($X_0, \mu_0$) is locally topologically equivalent to

Further, using a coordinate transformation and time reparametrization

system Eq (3.2) becomes

where

Since

According to the theory on the BT bifurcation in ^[42] and the preceding analysis, we have the following.

Theorem 3.1. Let $\lambda_1 = V_{Leak} + 63.413152$ and $\lambda_2 = g_{Na} - 2476.1036$. If bifurcation parameters $(V_{Leak}, g_{Na})$ vary around $(-63.413152, 2476.1036)$, then the system (6), at the Bogdanov-Takens bifurcation point BT, is locally topologically equivalent to the following system:

System (8) has the following local representation of bifurcation curves in a small neighborhood of the origin:

(1) There is a fold bifurcation curve

(2) There is a nondegenerate Hopf-bifurcation curve

(3) There is a saddle homoclinic-bifurcation curve

In Section 3.3, we conduct a bifurcation analysis of the model, with a specific focus on the BT bifurcation of the system. Through calculations on the simplified model, we obtain an in-depth analysis of the opening and breaking behavior associated with the BT bifurcation. Furthermore, we clarify the topological structure near the bifurcation point, particularly by calculating the homoclinic orbit curve and resolving issues related to its computation.

4.   Discussion and conclusions

In this paper, the kinetic properties of a class of dopamine neurons were investigated and abundant results were obtained. Due to the complexity of neurons, neuronal models are often characterized by high latitude and high complexity. We mainly focused on the dynamical properties of high-dimensional systems squeezing a reduced dimensional system. It is very difficult to study the dynamics of the model directly, so we first use the correlation of variables to downsize the original system in two steps and obtain a simplified system in three dimensions.

The first step is to consider the dimensionality reduction using the membrane potential. By drawing the phase diagram of each gating variable with respect to the membrane potential and substituting the appropriate data fitting, a simplified system in eight dimensions is obtained and the system retains the characteristic properties of the original system. The second step of dimensionality reduction is carried out with the help of the gating variables of the potassium channel $n$, and the calcium concentration $[Ca]$ to obtain the simplified model in three dimensions. It can be seen that the chaotic regions are generated due to the PD bifurcation, and the SN bifurcation on the LPC. Then, we study the equilibrium bifurcation of the system on the basis of a simplified model. The CP bifurcation, BT bifurcation, GH bifurcation, and other residual co-dimension-2 bifurcations of the simplified system are obtained through computation, the topology near the BT bifurcation is analyzed in detail, and the specific structure of each bifurcation on the central manifold near the neighborhood of this point is obtained.

To sum up, we give a method to study a high-dimensional nonlinear neuron model, and by comparing the action potentials generated by the model before and after simplification (the model after the second step of dimensionality reduction loses the firing characteristics of the original model), it is found that the chaotic regions are generated by the cascade formation of PD bifurcation generated by PD bifurcation (Figure 5, $V_{Leak} = [-60, -54]\; mV$), and the chaotic regions generated by SN bifurcation on the LPC ($V_{Leak} = [-56, -54]\; mV$). The inter-spike interval of the membrane potential of the system is studied to characterize the process of chaos generated by the PD cascade at $V_{Leak} = 58\; mV$. However, there is no further comparison of these differences in the generation of chaos, which can be further investigated. In addition, although the article gives the reason for the generation of chaos by combining co-dimension-2 bifurcation, it is difficult to analyze the dynamics of the original high-dimension neuron model, and the generality of the model simplification needs to be further proved.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

This work is supported by National Natural Science Foundation of China (No. 12072236), Henan Province Postgraduate Education Reform and Quality Improvement Project (No. YJS2023SZ15), and The key program of Henan province higher education (No. 24A110005)

Conflict of interest

The authors declare there is no conflicts of interest.