Reduced and bifurcation analysis of intrinsically bursting neuron model

1.   Introduction

The human nervous system represents one of the most intricate and complex systems in the body, forming the foundation of our perception, cognition, learning, and action. The fundamental electrical activity of neurons forms the basis of nervous system function. The interactions between neurons can form neural networks that support important functions including sensation, movement, learning, and memory^[1,2]. Neuron dysfunction can lead to a range of neurological disorders, including severe depression^[3,4], and Parkinson's disease^[5,6]. Consequently, extensive research has been conducted in the field of neuroscience investigating the electrical activity of fundamental neurons and the relationship between neural networks and behavior. These investigations are significant in comprehending the functional mechanisms of the nervous system and in the treatment of related neurological disorders.

Neurons, as the fundamental units of the nervous system, are key to understanding the physiology and functional mechanisms of the nervous system. Communication and information exchange between neurons occur through a series of electrical and chemical signals. Neuronal firing refers to the process of membrane potential changes. In order to investigate neuronal properties in greater depth, mathematical models of the neuron's electro-chemical activity are typically used. The dynamic properties of neurons can be described by a set of differential equations. These mathematical models provide the theoretical framework for investigating neuronal interactions and developing computer simulation programs to mimic the behavior of single neurons and neuronal networks^[7,8]. By exploring neuronal models and related computer programs, we can gain a deeper understanding of the intrinsic mechanisms underlying neuronal activity, assisting in the development of theoretical foundations for the entire nervous system.

The electrical activity of neurons is complex and diverse, exhibiting various forms of firing such as spiking, bursting, and mixed-mode oscillations^{[9,10,11,12,13]}. These rhythms generated by such electrical activity are also highly diversified, ranging from chaos to period-doubling bifurcations, and period adding bifurcations^[9,14]. The mechanisms underlying the production and evolution of these rhythms involve the interaction among neurons and the characteristics of dynamical systems. Thus, the study of bifurcation phenomena is essential to understand various physiological functions displayed by neurons. Bifurcation analysis not only enables the comprehension of the electrical activity of neurons, but also provides a new perspective and useful tools for investigating complex issues related to other neural systems. For instance, in the modeling of neurological diseases like epilepsy through the construction of neural models and networks, bifurcation dynamics can be employed to analyze the stability of the models and to investigate the impact of parameter changes on the progression and treatment of neurological disorders.

The dynamical models of neurons are often characterized by high-dimensional and highly nonlinear features, while bifurcation theory is commonly used for low-dimensional systems, such as two-dimensional or three-dimensional dynamical systems. Therefore, when studying high-dimensional neuronal systems, the problem of dimensionality reduction inevitably arises. Additionally, many variables in neuronal models are unobservable, rendering parameter estimation in neuronal models infeasible. Moreover, because the complex dynamics and processes of neuronal activity involve numerous details and factors, it is challenging to analyze and characterize them in a simple manner. These factors contribute to the high dimensionality of neuronal models and necessitate dimensionality reduction techniques for better understanding and analysis. Common dimensionality reduction methods are: the dimension reduction method based on the central manifold theory ^[15,16], Lyapunov-Schmidt (LS) method ^[17], Nonlinear Galerkin method ^[18,19], Proper Orthogonal Decomposition (POD) method ^[20,21], fast and slow dynamics analysis method and geometric singular perturbation theory^{[10,11,13,14]}.

The focus of this study is the intrinsic bursting neuron model, also referred to as the Pospischil model after its creator^[22]. Intrinsic bursting neurons are a prevalent type of nerve cell that can generate bursts of activity in response to depolarizing stimuli, with accompanying frequency adjustments. This firing pattern has been observed in both in vivo and ex vivo experiments in the primary sensory cortex. Experiments in guinea pig somatosensory cortex have recorded these patterns, as have studies of the core visual cortex in cats^[23,24]. Similar patterns have been observed in neocortical neurons as well, where bursting due to frequency adjustment after spike under depolarization has been noted. Mathematical modeling can aid in the understanding of the intrinsic bursting neurons. In this paper, we investigate the dynamic properties of a seven-dimensional Pospischil model, which is first reduced in dimensionality to facilitate analysis. Through correlation analysis between channel gating variables, the Pospischil model is reduced to a three-dimensional model, which is subjected to bifurcation analysis to explore the bursting mechanism of the model.

The rest of this paper is organized as follows. Section 2 presents the model and method. The results are provided in Section 3. Section 4 presents the conclusion and discussion.

2.   Model and methodology

High-dimensional Pospischil neurons studied in this paper are based on the mathematical model of ion channel, which is based on the minimum set of voltage-dependent conductance. In order to enable it to generate bursting behavior, based on the RS model, according to the experiment ^[25] the $L$-type calcium current is added to the model. The Pospischil model is a single-compartment neuron model whose membrane potential and Gates equations are

where $V$ is the membrane potential, $C_m$ is the membrane capacitance, $I$ denotes the external stimulus, $g_{leak}$ denotes the leak conductance, $V_{leak}$ is the leak reversal potential, and these parameters depend on the input resistance. $I_{Na}$ and $I_{Kd}$ are the sodium ion current and potassium current which are responsible for generating the action potential. $I_M$ and $I_L$ are the slow non-inactivation of potassium current and the high threshold calcium ion current, respectively ^[26,27,28], which is expressed as:

where $\alpha_x(V)$ and $\beta_x(V), (x = m, h, n, q, r)$ indicate the rates of gating of various ion channels from inactivation to activation and from activation to deactivation, respectively, and they are related to the variable $V$. $m$ and $h$ denote the activation gated and inactivated gated variables of the sodium ion channel, respectively. $n$ denotes the gated variable of the potassium ion. $q$ and $r$ denote the activation gated and inactivated gated variables of the calcium ion. The sodium ion gated conversion rate is expressed as:

where $V_T$ is the parameter that adjusts the firing threshold. Steady state and the time constant of $x$ can be written as:

The system (2.1)–(2.7) can be used to describe the fire activity of the intrinsic bursting neuron, and it can also produce the corresponding intrinsic firing characteristics, where $g_{Na}, \; g_{Kd}, \; g_M$ and $g_L$ are the maximal conductance for the various types of channels, $V_{Na}, \; V_K$ and $V_{Ca}$ are the reversal potential for the various types of channels. The specific parameter values of the model are shown in the Table 1, The temperature corresponding to all kinetic parameters is 36 $^{\circ}$C. The default values of the parameters are taken from the values in the Table 1.

In this paper, the model is numerically calculated using Python software and fourth-order Runge-Kutta methond, and the step is 0.001 ms. NEURON is used to process the correlation between gated variables^[29], Matcont^[30] and XPPAUT^[31] are used to the bifurcation analysis.

3.   Results

3.1.   Reducing the model

3.1.1.   Projection reduction method

In order to be able to use fewer ion channels to form the Pospischil model, the original 7-dimensional model needs to be dimensionally reduced. The idea of dimensionality reduction is from the Hodgkin-Huxley model to the simplified method of the FitzHugh-Nagumo model^[32]. By NEURON, the data of various gated variables at the same time can be obtained. By processing the data, a definite functional relationship with each other can be found. Using these relationships, one variable can be represented by additional variables, thereby the dimensionality of the system of equations can be reduced. The operation method is reproduced below. Consider $n$ dimensional differential equations:

Suppose there is a functional relationship between $x_i$ and $x_j$, $1\leq i, j\leq n, \; i\neq j$, i.e.,

Substitute (3.2) into (3.1), then the original model becomes an algebraic differential equation:

This means using (3.2) to project the system (3.1) to $(n-1)$ dimension space spanned by $\{x_1, x_2, \cdots, x_{j-1}, x_{j+1}, \cdots, x_n\}$. So the reducing the model can be used to find the relationship between different variables.

3.1.2.   Reducing Pospischil model

In the neuron model, the relationship between the gated variables reflects the regulatory mechanism of the corresponding ion channel. Therefore, the model can be simplified by quantifying the regulatory relationship between the ion channels. Figure 1b gives the relationship between the potassium channel gating variable $n$ and the sodium ion channel inactivation gating variable $h$. Combined with the regulation mechanism of potassium ion and sodium ion gating variables, the seven-dimensional Pospischil model is reduced to a six-dimensional model by using the projection dimensionality reduction method, and the similarities and differences between the reduced-dimensional model and the original model are compared.

The seven-dimensional Pospischil model (2.1)–(2.7) is

Let's consider the gating variable $n$ in the Eqs (2.1)–(2.7), by means of the control mechanism between different gating variables, making it a six-dimensional model. Then compare the similarities and differences between the firing generated by the model after dimension reduction and the the original model.

First, we consider that the regulation mechanism between the sodium ion inactivation gated variable $h$ and the potassium ion gated variable $n$ is linear, so that we can get:

The diagram of Eq (3.3) is the black straight line in Figure 1b). A simplified six dimensional model can be obtained by substituting Eq (3.3) into the Pospischil model (2.1)–(2.7). When $g_L = 0.365\; \rm{mS/cm^2}$, The firing of the new model is shown in Figure 1c).

Then a more complex regulatory mechanism is considered. If the sodium ion inactivation gating variable $h$ and the potassium ion gating variable $n$ change as follows:

The diagram of Eq (3.4) is the purple curve in Figure 1b). The adjusted fitting coefficient is 0.9666. A simplified six dimensional model can be obtained by Substituting Eq (3.4) into the Pospischil model (2.1)–(2.7). When $g_L = 0.165\; \rm{mS/cm^2}$, The firing of the new model is shown in Figure 1d).

Next the regulatory mechanism in an exponential form is considered:

where $a = 0.8956, \; b = -2.487, \; c = -2.174\times 10^{-14}, \; d = 28.57$. The diagram of Eq (3.5) is the light blue curve in Figure 1b). When $g_L = 0.175\; \rm{mS/cm^2}$, The firing of the new model is shown in Figure 1e). From the data of $(h, n)$ in Figure 1b), the regulatory relationship between $h$ and $n$ is approximated as an ellipse. So just try to use elliptic functions to control the gated variables. The mechanism is characterized. Since the ellipse corresponds to a multi-valued function, the explicit relationship between $n$ and $h$ cannot be directly given, and the relationship between them is calculated by means of its first derivative. Suppose the variables $h$ and $n$ satisfy the function

Take the derivative of both ends of the Eq (3.6) respect to $t$, we can get

By mean of (2.3), (2.4), (3.6) and (3.7) explicit relationship between $h$ and $n$ can be obtained:

Take the derivative of both ends of Eq (3.8) respect to $t$, we obtain:

Substituting Eqs (2.3) and (2.4) into (3.9), we can get

The final relationship between $n$ and $h$ can be obtained by simultaneous Eqs (3.8) and (3.10):

with $B = 1, C = 1, D = -2, E = -2, F = 1.$ When $g_L = 0.095\; \rm{mS/cm^2}$, the time course of reduced model can be shown in Figure 1f), and the fitting curve is the green curve in Figure 1b).

The present study investigated the regulation of bursting activity in seven-dimensional neuronal models through a multi-valued control mechanism. Our analysis revealed that the essential features of basic bursting neurons, such as the generation and frequency adjustment of bursts during depolarization, were retained regardless of the regulatory mechanism employed. Moreover, by tuning the parameter $g_L$, the form of the bursting, including its start time, duration, and number of spikes, could be consistently maintained through all models. However, the shape of a spike was found to differ from that produced by the original system, owing to the use of multi-valued functions to regulate. Specifically, we found that the Pospischil model can be effectively simplified by reducing its dimensionality without increasing its nonlinearity, while retaining all intrinsic bursting characteristics, through the use of linear control between two gated variables.

In conclusion, Figure 1 illustrates the period-9 bursting produced by the original seven-dimensional model, as well as the bursting profiles produced by simplified models under various regulatory mechanisms. Our work highlights the utility of multi-valued control mechanisms in regulating bursting activity in neuronal models, while maintaining essential features of such neurons. Specifically, our findings suggest that simplification of neuronal models can be achieved by reducing their dimensionality through linear regulatory mechanisms, without compromising on their intrinsic bursting features. These insights could inform the development of simplified neuronal models with applications in the study of neuronal dynamics and in neural engineering.

Step 1: Collecting information on the associated ion channel gated variables. The NEURON software is used to process the model and record the time course of gated variables of the ion channel. The results are shown in Figure 2. Figure 2a)–c) gives the relationship between the gated variables $h$ and $q$, $m$ and its steady state $m_{\infty}$ and $p$ and $r$, respectively.

Step 2: Combined with the Figure 1b), the relationship of the relevant gated variables can be obtained:

where $c_0 = 0.85, \; c_1 = 0.95, \; c_2 = -0.3243$ and $c_3 = 0.5485$.

Step 3: Replace $n$ and $q$ in the original seven-dimensional Eqs (2.1)–(2.7) with $h$ by the projection reduction method, and $m$ with $m_{\infty}$. And $r$ is replaced with $p$, so the system (2.1)–(2.7) is reduced to a three-dimensional model:

The specific ion channel models are:

where $c_0 = 0.85, \; c_1 = 0.95, \; c_2 = -0.3243$ and $c_3 = 0.5485$. The bursting generated by system (3.11)–(3.13) is shown in Figure 2d). It can be seen from the figure that the bursting generated by the reduced system still retains the basic characteristics of the intrinsic bursting neurons, and the shape of the burst is basically the same as that of the seven-dimensional original model. Therefore, the reduced three-dimensional model is studied below.

3.2.   Codimension 2 Bifurcation of reduced Pospischil model

The present study focuses on the analysis of dynamical properties of the Pospischil model (3.11)–(3.13) after dimensionality reduction, which facilitates a better understanding of its dynamics. To delve further into the dynamics of the reduced model, we introduce the consideration of the slow non-inactivation potassium current, which when added to the model, enables the adjustment of the firing frequency. In this context, the maximum conductance of the ion channel $(g_M)$ is selected as the primary parameter of interest in studying the dynamics of the three-dimensional model. Additionally, we examine the impact of varying external stimuli on the firing pattern of the model by considering $I$ as a second parameter.

Building on the foundation of the equilibrium point single-parameter bifurcation diagram, we extend our analysis to include a bi-parameter bifurcation diagram of equilibrium points for the reduced system in the parameter plane $(g_M, \; I)$ as shown in Figure 3.

It's shown in Figure 3, a bi-parameter bifurcation diagram of the equilibria of the reduced model. The abscissa is $g_M$ and the ordinate is $I$. On the $(g_M, \; I)$ plane, the red curve represents the saddle node bifurcation. They are denoted as $SN_1$ and $SN_2$ respectively. The blue curve represents the Hopf bifurcation, denoted as $H_1$ and $H_2$ respectively. There is a neutral saddle on the curve, and the corresponding parameter coordinates are $(0.36467, \; 2.54933)$. Two Bogdanov-Takens bifurcations appear at the intersection of the saddle-node bifurcation curve $SN_1$ and the Hopf curve $H_1$, denoted as $BT_1$ and $BT_2$, respectively, and the corresponding parameter coordinates are $BT_1 \; (0.002628, 0.043215)$ and $BT_2 (-4.930709, -0.285314)$. Following consideration of the physiological significance of the parameter $g_M$, we study the dynamical behavior of the system at $BT_1$. Figure 3b) is a partial enlargement of Figure 3a). According to the bifurcation theory ^[16], in addition to the saddle note curve and the Hopf curve at the Bogdanov-Takens bifurcation, there will be a homoclinic bifurcation curve, but because the distance between the saddle-note curve $SN_1$ and the Hopf curve $H_1$ is very small, it is difficult to describe the homoclinic bifurcation curve with the accuracy of the existing computer. Therefore, the bifurcation behavior of the reduced system at the Bogdanov-Takens bifurcation $BT_1$ is analysed on its central manifold^[33].

By calculation, $BT_1$ can be gotten in the Figure 3, where $(g_M, I)^T = (0.002628, 0.043215)^T \triangleq \mu_0$, and the corresponding equilibrium point coordinates are $(V, h, p)^T = (-54.898367, 0.995835, 0.120274)^T \triangleq X_0$, rewrite the system (3.11)–(3.13) as:

where $X = (V, h, p)^T, \mu = (g_M, I)^T$, and

The Taylor expansion of $F(X, \mu)$ around $(X_0, \mu_0)$ is

where

The eigenvalues of matrix $A$ are 0, 0, and $-0.3291554654$.

Let $P = (p_1, p_2, p_0)$ and $P$ is an invertible matrix. We change matrix $A$ into the Jordan normal form. we can obtain

with

So

Let $P^{-1} = (q_1, q_2, q_0^T)^T$. Then we have

By a straightforward calculation as described, we obtain

Note that $\lambda_1 = g_M-0.002628$ and $\lambda_2 = I-0.043215$. It gives

So the dynamics on the center manifold of system (3.11)–(3.13) is locally topologically equivalent to the following system at the Bogdanov-Tankens bifurcation $(X_0, \mu_0)$:

Furthermore, by virtue of the transformations

the system (3.15) can be rewritten to

where $\bar{\beta_1} = -2.432461968\times 10^{5}\lambda_1+57616.51458\lambda_2$ and $\bar{\beta_2} = -1059.090571\lambda_1-77.07800081\lambda_2$.

Applying the theory on the Bogdanov-Tankens bifurcation^[33], we get

Theorem 3.1. Let $\lambda_1 = g_M-0.002628$ and $\lambda_2 = I-0.043215$. If the values of parameters $(g_M, \; I)$ changes around $(0.002628, \; 0.043215)$, and other parameters in system (3.11)–(3.13) are indicated as in Section 2. Then system (3.11)–(3.13) is locally topologically equivalent to the following system at the Bogdanov-Tankens bifurcation ${\rm BT_1}$:

which has the following local bifurcation behaviors in a small neighborhood of the origin:

(i) There is a saddle-node bifurcation curve

(ii) There is an Andronov-Hopf bifurcation curve

(iii) There is a homoclinic bifurcation curve

This method can be used at the Bogdanov-Tankens bifurcation ${\rm BT_2}$.

4.   Discussion

This article presents a study on the dimensionality reduction of high-dimensional neuron models and their dynamical oscillation phenomena. Specifically, the Pospischil neuron model is analyzed to demonstrate the dimensionality reduction of the model along with bifurcation analysis and bursting mechanism after dimension reduction. It is important to note that the bursting neuron model is not limited to any particular neuron type, but rather describes the discharge characteristics common to a class of neurons found in various physiological experiments on the cerebral cortex. These neurons are also major in constructing thalamic-cortical circuits. In modeling such neurons, they are powered by a seven-dimensional system description taking into account various ion channels. However, since it is difficult to theoretically study high-dimensional dynamic systems, it is necessary to streamline the model.

In this study, we begin with a seven-dimensional Pospischil model and apply a method of dimensionality reduction that differs from previous methods. Specifically, we utilize a method that examines the correlation between discrete gated variables in the model. By comparing various ion channels gated variables, we identified a regulatory mechanism that allows us to simplify the model using a simple linear relationship. As a result, the original model was reduced to a three-dimensional model by considering $h$ and $n$, $h$ and $q$, $p$ and $r$, and $m_{\infty}$ and $m$ as a single variable. Our reduced model preserves the intrinsic bursting characteristics of the original model, which includes burst under depolarization stimulation, discharge, and frequency adjustment. Therefore, we are able to make a simplified single parameter adjustment (changing the value of $g_L$) to ensure that the reduced model produces a substantially uniform bursting that maintains a considerable degree of consistency regardless of the amplitude of the firing, the frequency of the bursting, and the duration. Although the three-dimensional model cannot fully inherit the complex dynamic phenomena in the original model, we find that such dimensionality reduction methods have certain universality and provide an opportunity for studying high-dimensional neuron models.

The present study investigates the bifurcation of equilibrium points for the reduced model, focusing on bi-parameter bifurcation. As the calcium ion channel and external stimulus act as the primary factors in generating bursting, the maximum conductance of slow non-inactivated potassium channel ($g_M$) and external stimulus ($I$) are explored as bifurcation parameters. Moving on to bi-parameter bifurcation analysis, we investigate the system's generation of Bogdanov-Takens bifurcation. Theoretical analysis reveals that the system is generated in close proximity to the Bogdanov-Takens bifurcation, and conditions for saddle bifurcation, Hopf bifurcation and homoclinic bifurcation are derived. Numerical simulation is used to verify the correctness of the theory and analysis.

In the study of the model, not only the characteristics of the intrinsic bursting neurons but also the bifurcation characteristics of other neuronal types can be observed. However, it is imperative to note that dimensionality reduction techniques could suffer from the issue of information loss. For example, the ISIs of sequence of the membrane potential in the two systems is basically the same, and its main features–double period bifurcation, add-period bifurcation and chaotic area–are present in both systems, indicating the reduced model can better reflect the essential features of the original model. However, a little difference can be observed in the Figure 4. First, the period-doubling bifurcation of the original model appears at $g_M = 0.0338\; {\rm mS/cm^2}$ Nearby, the chaotic area appears in $g_M = 0.03425\; {\rm mS/cm^2}$ Recently. Second, the original model forms a chaotic region at the same time when the first period of bifurcation is generated, but the chaotic region is simplified in the reduced model. This may be explained by the complexity of the high-dimensional model itself, so this property is difficult to observe after dimension reduction. Looking for the hidden chaotic attractors is the future work^[34,35,36].

The autonomous neuronal model is considered in this paper, but for the non-autonomous systems, the incremental integral transformation is the appropriate method^[37,38,39]. Therefore, reducing the dimensionality of neuronal dynamics models can be challenging. This challenge arises since it is necessary to maintain a detailed description of specific behaviors while attempting to simplify the model for better operability and comprehensibility. Thus, discovering a reasonable approach for effectively reducing the dimensionality of neuronal dynamics models for improved analysis and comprehension of their behavior is an essential facet of studying neurons. Nonetheless, the precise influential mechanisms of the reduction practices require further scholarly investigation.

Use of AI tools declaration

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

Acknowledgments

This work was supported by the National Science Fund for Distinguished Young Scholars (No.12101198), and Research Project of Henan province postdoctoral (No. 19030095).

Conflict of interest

All authors declare no conflicts of interest in this paper.