Approximate solution of nonlinear fuzzy Fredholm integral equations using bivariate Bernstein polynomials with error estimation

1.   Introduction

Individual neurons are the basic units of the biological nervous system, they process information by working together in neuronal circuits with specific synaptic connectivity. In response to an applied current, the signal can be propagated from one neuron to the next neuron by generating electrical impulses or spikes^[1]. Hodgkin and Huxley (HH) firstly proposed the mathematical description of the neural model based on a large number of biological experiments^[2], which built a new field of neuronal electrophysiology and laid the foundation for theoretical neuroscience. However, the HH model is a complex nonlinear system, and is difficult to be analyzed theoretically. Subsequently, researchers have successively proposed many simplified neuron models such as the Fitzhugh-Nagumo (FN) model ^[3], Morris-Lecar (ML) model^[4], Hindmarsh-Rose (HR) model^[5], and the modified forms of them^[6,7,8], to explore the physiological mechanisms of neurons.

Previous studies have shown that many neurological dynamic diseases such as Epilepsy, Parkinson's disease, Alzheimer's disease and Schizophrenia may be induced by bifurcations caused by the changes of the regulatory parameters in the neural system, which leads to the disorders in the firing patterns of nervous systems^[9,10,11,12]. In addition, the changes of firing patterns in neuron systems may cause neurons to produce a kind of plasticity, which finally affects learning, memory and other neural activities that rely on neuron-plasticity^[13]. Therefore, a method is expected to be proposed to control the bifurcation mechanism of the neuron model, which can change the firing patterns from abnormal to normal under certain conditions.

Bifurcation control has attracted increasing attention due to its potential applications in many fields, such as engineering, biomedicine, and meteorology systems etc^[14]. The main idea of bifurcation control is to change the inherent bifurcation characteristics of the original system by designing the controller without changing the parameters, thereby obtaining some desirable bifurcation behaviors. This includes delaying or advancing the onset of the bifurcation points, modifying the shape or type of bifurcations, and so on. Many bifurcation control methods have been provided in previous studies^{[15,16,17,18,19]}. Xie et al. ^[10] studied the dynamic behaviors of Hopf bifurcation under the control of the washout filter in HH model. The results show that the linear term of the controller can control the position of the Hopf bifurcation point, and the cubic term can change the criticality of the Hopf bifurcation point. The washout filter dynamic feedback controller is introduced into the two-dimensional HR model to change the firing pattern from type I to type II by delaying the location of its inherent Hopf bifurcation^[12]. Huang et al. ^[19] studied the Hopf bifurcation control of the ML model with type I under the control of washout filter-aided dynamic feedback controller. However, the firing pattern cannot be changed from type I to type II by using the washout filter-aided dynamic feedback controller in ML model. Therefore, a new controller should be introduced to address this issue. In addition, the analytical expression for determining the stability index of Hopf bifurcation criticality has not been obtained. Moreover, none of these studies involved the bifurcation control of the ML model with Hc bifurcation type.

Motivated by above analysis, in the present paper, we introduce a new controller proposed by Nguyen et al.^[20] into the ML model to study the influence of the controller on the ML model with Hc bifurcation type. In addition, we give the analytical expression of stability index of the system to determine the criticality of Hopf bifurcation.

The structure of the paper is organized as follows. In Sect. 2, we reviewed the bifurcation dynamics and the conditions for the emergence of Hopf bifurcations of the ML model. In Sect. 3, the dynamic state-feedback control law was applied to the ML neuron model with Hc bifurcation type. Also, we derived the analytical expression of the bifurcation stability index for the controlled ML model with Hc bifurcation type. Conclusions are given in Sect. 4.

2.   Reviews of the ML model

2.1.   ML model

The ML neuron model is a simplified form of HH model and its mathematical description can be expressed by the following nonlinear differential equations:

with

where $V$ is the membrane potential, and $N$ is the recovery variable. $I$ is the current stimulation and is treated as the bifurcation parameter. $C$ is the membrane capacitance. $V_{i} \ (i = K, Ca, L)$ represent the reversal potential, and $g_{i} \ (i = K, Ca, L)$ represent the maximum conductance. $V_{i} \ (i = 1, 2, 3, 4)$ are constant potentials, and $\phi$ represents the temperature factor. Here, we named system (1) as the original system.

2.2.   The bifurcation dynamics of the ML model

There are three types of bifurcations in ML model by choosing different parameter sets, including Hopf bifurcation, saddle node bifurcation on invariant cycle (SNIC) and homoclinic (Hc) bifurcation (see Table 1). In this paper, we mainly focus on the Hopf bifurcation control of the ML model with Hc bifurcation type. The bifurcation diagram of the ML model with Hc bifurcation type is shown in Figure 1. From Figure 1 we can see that as the bifurcation parameter $I$ increases, the ML model with Hc bifurcation type is divided into six types of dynamic regions. In region I (when $I < -9.949$), the system has only one stable equilibrium point in the lower branch; In region II (when $-9.949 < I < 34.94$), the system contains three equilibrium points, the stable equilibrium point in the upper branch, the unstable equilibrium point in the middle branch and the stable equilibrium point in the lower branch; In region III (when $34.94 < I < 36.32$), the system contains a stable limit cycle and three equilibrium points; In region IV (when $36.32 < I < 39.96$), the system contains a stable limit cycle, an unstable limit cycle and three equilibrium points; In region V(when $39.96 < I < 40.56$), the system contains a stable limit cycle, an unstable limit cycle and a stable equilibrium point; In region VI(when $I > 40.56$), the system contains a stable equilibrium point. The corresponding phase diagrams of six typical dynamical behaviors are shown in Figure 2. Compared it with ^[19], the ML model with Hc bifurcation type has more dynamical behaviors than that with SNIC bifurcation type. In particular, Hc bifurcation occurs in the system when $I = 34.94$, and the corresponding phase diagram is shown in Figure 3.

2.3.   Two types of firing patterns in ML model

Some researchers have proposed that the firing patterns can be determined by the types of bifurcation^[21,22]. Specifically, type I firing pattern occurs when the system undergoes SNIC bifurcation or Hc bifurcation, while type II firing pattern occurs when the system undergoes Hopf bifurcation. For type I, the system can fire at any low frequency; for type II, the system fires from a positive minimum frequency. Figure 4 shows the firing patterns of the ML model under the three types of bifurcations. From Figure 4 we can see that the ML model with Hopf bifurcation type shows the type II firing pattern (see Figure 4b), and the corresponding time evolutions are shown in Figure 4a; the ML model with SNIC bifurcation or Hc bifurcation type both show the type I firing pattern (see Figure 4d and 4f), and the corresponding time evolutions are shown in Figure 4c and 4e, respectively.

2.4.   Conditions for the emergence of Hopf bifurcation

Consider the following n-dimensional nonlinear system:

where $x\in R^{n}$ represents the state variables, $\mu \in R$ represents the bifurcation parameter. Assume that the system has a fixed point $x_{0}$ when the bifurcation parameter $\mu = \mu_{0}$ and the corresponding eigenvalue $\lambda_{\mu_{0}} = \alpha_{\mu_{0}}+i\beta_{\mu_{0}}$. The Hopf bifurcation occurs at $\mu = \mu_{0}$ when the system satisfies the following conditions: $\beta_{\mu_{0}}\neq0$, $\alpha_{\mu_{0}} = 0$, and $\frac{d}{d\mu}(\alpha_{\mu})\mid\mu = \mu_{0} > 0$. However, the eigenvalues of the system cannot be obtained easily when the system is a high-dimensional system. Therefore, the following Routh-Hurwitz stability criterion can be used to determine whether the Hopf bifurcation occurs^[23]. Assume that the corresponding characteristic polynomial of system (3) is as follows:

Then, we obtain

where $p_{i}(\mu) = 0$, if $i < 0$ or $i > n$. Then, we define

Based on the Routh-Hurwitz stability criterion, the Hopf bifurcation occurs at $\mu = \mu_{0}$ when the system satisfies the following conditions:

(I) Eigenvalue crossing condition:

(II) Transversality condition:

3.   Hopf bifurcation control of the ML model

3.1.   The dynamic state-feedback controller

In this section, we provide a dynamic state-feedback controller to control the Hopf bifurcation of the original system. Consider the following nonlinear system:

where $y\in R^{m} (1 \leq m \leq n)$ represents the state variables of the controller, $u(x, y)$ represents the feedback control and $g(x, y)$ is a smooth function. In this paper, we utilize the following feedback control law proposed by Nguyen and Hong^[20].

where $x_{i}^{0} (i = 1, 2, ..., m)$ are the equilibrium points of the system at the Hopf bifurcation point, $k_{1i}$ and $k_{3i}$ represent the linear term control gain and the cubic term control gain, respectively. $l_{i}$ are constant parameters. Then we obtain the following controlled system:

where

According to Eq (10), we have $u_{i} = 0$ when $\dot{y} = u_{i} = 0$. So, if $x^{0}$ is the equilibrium point of the original system, $(x^{0}, y^{0})$ is the equilibrium point of the controlled system, where $y^{0} = (y^{01}, y^{02}, \cdots, y^{0m})$ and then $y^{0i} = (k_{1i} x_{0i}+k_{3i}(x_{0i}-x_{01i})3)/l_{i} (i = 1, 2, \cdots, m)$. That is, the equilibrium points of the original system are the equilibrium points of the controlled model. The reason for using state feedback control to conduct the Hopf bifurcation control is that the equilibrium structure of the original system and the controlled system remains unchanged during the control process. In the present paper, we add the controller to the membrane potential V, and then we obtain the following controlled system:

where $V^{0}$ is the equilibrium membrane potential of the original system at the Hopf bifurcation point.

3.2.   Hopf bifurcation control of the ML model with Hc bifurcation type

The Hopf bifurcation control of neuron models has been widely studied, including HH model, HR model, ML model with Hopf bifurcation type or SNIC bifurcation type^{[10,11,12,19,20]}. However, the Hopf bifurcation control of the ML model with Hc bifurcation type has not been studied. In addition, the controllers used in previous researches only studied the system with single Hopf bifurcation point. Therefore, in the present paper, we utilize the dynamic state-feedback control proposed by Nguyen et al. to make the Hopf bifurcation point change from one to two. The parameter values of the model are $\phi = 0.23, g_{Ca} = 4, V_{3} = 12, V_{4} = 17.4, V_{Ca} = 120, V_{K} = -84, V_{L} = -60, g_{K} = 8, g_{L} = 2, V_{1} = -1.2, V_{2} = 18, C = 20, l = 0.1$, and $I$ is treated as the bifurcation parameter.

3.2.1.   The transition from Hc bifurcation type to SNIC bifurcation type

The Hopf bifurcation point is located at $I = 36.32$ in the original ML model with Hc bifurcation type (Figure 1). To change the model from the Hc bifurcation type to the SNIC bifurcation type, the Hopf bifurcation point needs to be moved forward. Here, we change the Hopf bifurcation point from $I = 36.32$ to $I_{0} = 70$. The equilibrium point of the original system is $(V^{0} = 6.7697, N^{0} = 0.35407)$. Therefore, the equilibrium point of the controlled system is $(V^{0} = 6.7697, N^{0} = 0.35407, y^{0} = 67.697k_{1})$, and then we obtain the Jacobian matrix of the controlled system at the equilibrium point:

Then, we obtain the characteristic polynomial of the Jacobian matrix $J(I_{0})$:

where $p_{0} = 1, p_{1} = 0.1956-k_{1}, p_{2} = 0.1992-0.2326k_{1}, p_{3} = 0.01896$. Substituting $p_{i}\ (i = 0, 1, \cdots, m)$ into the crossing condition (I), we have

Solving Eq (16), we have $k_{1} = 0.0893$. Then, we have:

which satisfy the transversal condition (II).

The bifurcation diagram of the controlled system when $k_{1} = 0.0893, k_{3} = 0$ is shown in Figure 5a. As expected, the Hopf bifurcation point of the ML model with Hc bifurcation type moves forward from $I = 36.32$ to $I = 70$. Figure 5b shows the two-parameter bifurcation diagram of bifurcation parameter $I$ with respect to $k_{1}$, which further verifies that Hopf bifurcation occurs at $I = 70$ when $k_{1} = 0.0893$. Now, the model shows a change from the Hc bifurcation type to the SNIC bifurcation type (see Figure 5a), but does not change the firing pattern of the system (see Figure 5c and 5d).

The results of this section show that the ML model can be transformed from the Hc bifurcation type to the SNIC bifurcation type by introducing a dynamic state-feedback control. In addition, there is only one Hopf bifurcation point in the system when $k_{1}\neq0, k_{3} = 0$ (Figure 5a), that is, $k_{1}$ can change the location of the Hopf bifurcation point of the ML neuron model with Hc bifurcation type, resulting in the change of bifurcation type from Hc type to SNIC type.

3.2.2.   The transition from Hc bifurcation type to Hopf bifurcation type

Previous study has investigated the transition of ML model from SNIC bifurcation type to Hopf bifurcation type^[18]. Here, we will study the transition of the ML model from Hc bifurcation type to Hopf bifurcation. First, the inherent Hopf bifurcation point in the Hc bifurcation type needs to be moved forward, and then a new Hopf bifurcation point is constructed in front of the moved inherent Hopf bifurcation point. Here, we expect to change the inherent Hopf bifurcation point from $I = 36.32$ to $I_{1} = 200$, and then a new Hopf bifurcation point is constructed at $I_{2} = 60$.

When $I_{1} = 200$, the equilibrium point of the controlled system is $(V^{0} = 12.94, N^{0} = 0.52697, y^{0} = 129.4k_{1})$, and then we obtain the Jacobian matrix of the controlled system at the equilibrium point:

Then, we obtain the characteristic polynomial of the Jacobian matrix $J(I_{1})$:

where $p_{0} = 1, p_{1} = 0.4676-k_{1}, p_{2} = 0.3243-0.2301k_{1}, p_{3} = 0.028756$. Substituting $p_{i} \ (i = 0, 1, \cdots, m)$ into the crossing condition (I), we obtain:

Solving Eq (20), we have $k_{1} = 0.34966383$. Then, substituting $k_{1} = 0.34966383$ into the controlled system when $I_{2} = 60$, we obtain the equilibrium point of the controlled system $(V^{0} = 6.1368, N^{0} = 0.33762, y^{0} = 21.46+3148.761k_{3})$. The Jacobian matrix of the controlled system at the equilibrium point is as follows:

Then, we obtain the characteristic polynomial of the Jacobian matrix $J(I_{2})$:

where $p_{0} = 1, p_{1} = -0.1808-138.8506k_{3}, p_{2} = 0.1033-32.3938k_{3}, p_{3} = 0.0178$. Substituting $p_{i} \ (i = 0, 1, \cdots, m)$ into the crossing condition (I), we obtain:

Solving Eq (23), we have $k_{3} = -0.00205656$. Figure 6 shows the dynamic behaviors of the controlled system when $k_{1} = 0.34966383, k_{3} = -0.00205656$. From Figure 6a we can see that the Hopf bifurcation point of the ML neuron model with Hc bifurcation type is moved forward from $I = 36.32$ to $I_{1} = 200$, and a new Hopf bifurcation point is constructed at $I_{2} = 60$. Figure 6b shows the two-parameter bifurcation diagram of bifurcation parameter $I$ with respect to $k_{1}$, which further verifies that the model undergoes the Hopf bifurcation at $I_{1} = 200, I_{2} = 60$ when $k_{1} = 0.0893$. Figure 6c shows the time evolutions of the controlled system with different bifurcation parameter $I$. From Figure 6c we can see that the system fires from a positive minimum frequency, which is further verified in Figure 6d. Now, the model shows a change from Hc bifurcation type to Hopf bifurcation type (see Figure 6a), and the firing pattern changes from type I to type II (see Figure 6d).

The results in this section show that the ML model can be transformed from the Hc bifurcation type to the Hopf bifurcation type by introducing dynamic state-feedback control. In addition, the model performs type II firing pattern when $k_{1}\neq0, k_{3}\neq0$, that is, the combination of $k_{1}$ and $k_{3}$ can change the firing patterns of the ML neuron model with Hc bifurcation type.

3.2.3.   Theoretical analysis of the criticality of Hopf bifurcation

As mentioned above, the criticality of the Hopf bifurcation can be changed by a dynamic state-feedback controller. In this section, we derive the stability coefficient $\beta_{2}$ of the Hopf bifurcation by utilizing the Poincare Birkhoff (PB) normal form method. The Hopf bifurcation is subcritical when $\beta_{2} > 0$, resulting in an unstable limit cycle; the Hopf bifurcation is supercritical when $\beta_{2} < 0$, resulting in a stable limit cycle.

Consider the following controlled system:

where $x = (V, N, y)^{T}$, $f = (f_{1}, f_{2}, f_{3})^{T}$. The equilibrium point of the system at $\mu_{0}$: $x_{0}(\mu_{0}) = (V_{0}(\mu_{0}), N_{0}(\mu_{0}), y_{0}(\mu_{0}))^{T}$. Rewrite the controlled system as: $\dot x = f(x, \mu) = Jx+g(x)$, where $J$ is the Jocabian matrix of the system at the equilibrium point, $g(x)$ is the nonlinear term:

Then, we conduct the linear transformation of the controlled system: $x = Py+x_{0}(\mu_{0})$, where $P = [Re(v_{1}), Im(v_{1}), v_{3}]$, $Re(v_{1})$ and $Im(v_{1})$ are the real and imaginary parts of the eigenvector $v_{1}$ corresponding to eigenvalues $\lambda_{1}$, and $v_{3}$ is the eigenvector corresponding to the eigenvalue $\lambda_{3}$. Then, we have:

Substituting Eq (25) into Eq (26), we have:

where

$A = P^{-1}JP$ is the Jocabian matrix of the system, $F$ is the nonlinear term:

That is,

Based on the PB normal form theory, we obtain the expression of the bifurcation stability coefficient $\beta_{2}$:

where

Here, the vectors $w_{11}^{1}$ and $w_{20}^{1}$ satisfy the following linear equations:

where $\lambda_{3}$ is the real characteristic root of the Jacobian matrix, $\omega_{0}$ is the imaginary part of the pure virtual characteristic root of the Jacobian matrix, and $\omega_{0} > 0$. Solving Eqs (39) and (40) we obtain:

Then, we obtain the nonlinear terms of the controlled system:

Because $Ay$ is the first-order function with respect to $y$ and the second or more derivative of $y$ is required in the following calculations, the derivative of $y$ is omitted. We only calculate the following equation: $F = P^{-1}f(Py+x_{0}(\mu_{0}))$.

Next, we calculate the stability coefficient $\beta_{2}$ of the controlled system when $k_{1} = 0.3497$, $k_{3} = -0.0021$, $\mu_{0} = I_{1} = 200$. The equilibrium point of the controlled system is $(V_{0} = 12.94, N_{0} = 0.52697, y_{0} = 45.2465)$, and the corresponding Jacobian matrix is:

Then, we obtain the eigenvalues and eigenvectors: $\lambda_{1} = 0.4938i$, $\lambda_{2} = -0.4938i$, $\lambda_{3} = -0.1179$.

So, the matrix $P = [Re(v_{1}), Im(v_{1}), v_{3}]$ can be rewrite as:

Then, we obtain:

Substituting Eqs (48) and (49) into $F = P^{-1}f(Py+x_{0}(\mu_{0}))$, we have:

Finally, substituting the parameters into the stability coefficient formula we have $\beta_{2} = -9.9169\times10^{-4} < 0$, that is, the controlled system undergoes the supercritical Hopf bifurcation at $(V_{0}, N_{0}) = (12.94, 0.52697)$ when $I_{1} = 200$. The membrane potential curve and phase trajectory of the controlled system at the bifurcation point are shown in Figure 7a and 7b. From Figure 7b we can see that the phase trajectory of the controlled system converges to a stable limit cycle when $I_{1} = 200$. However, the phase trajectory of the original system converges to a stable equilibrium point when $I_{1} = 200$ (see Figure 7c and 7d). The similar results can be obtained when $I_{2} = 60$, so we omit them here.

4.   Conclusions

In this paper, we investigate the Hopf bifurcation control of the ML model with Hc bifurcation type by utilizing the method of dynamic state-feedback control. The results show that (1) the location of the Hopf bifurcation point can be changed by introducing the state-feedback control without changing the model parameters; (2) the control gain $k_{1}$ transforms the ML model from the Hc bifurcation type to the SNIC bifurcation type, but does not change the firing pattern of the model; (3) the cooperation of control gain $k_{1}$ and $k_{3}$ transforms the ML model from the Hc bifurcation type to the Hopf bifurcation type, and the firing pattern shows a change from Type I to Type II. In addition, we find that the criticality of the Hopf bifurcation of the system can be changed by introducing the control gain $k_{3}$, that is, the criticality of the Hopf bifurcation is changed from the subcritical Hopf bifurcation to the supercritical Hopf bifurcation. Many neurodegenerative diseases are caused by the change of firing patterns. Therefore, the results of this paper may have potential applications for the treatment of these diseases.

Compared to previous studies, the conclusions of this paper further enrich the research results. Some scholars have studied the Hopf bifurcation control based on the ML model with SNIC bifurcation type by introducing a washout filter-aided dynamic feedback controller^[19]. In their research, they only studied how the location and criticality of the Hopf bifurcation changed under different control conditions. However, they cannot perform the change of the Hopf bifurcation point from one to two. Therefore, it is impossible to study how to control the firing pattern of the system changes from type I to type II. Here, we utilize a new dynamic state feedback controller to control the Hopf bifurcation of the ML model with Hc bifurcation type. This controller not only reproduces the previous results, but also realizes the control of two Hopf bifurcation points. We further study the effects of the controller on firing patterns, providing potential application value for the treatment of some dynamic diseases. In addition, the stability index to determine the criticality of the Hopf bifurcation is given in our paper. The mutual verification of theoretical analysis and numerical simulation improves the integrity of this research.

Acknowledgments

This work was supported by the national Natural Science Foundation of China under Grant 11972115 and the Fundamental Research Funds for the Central Universities. The authors would like to thank the reviewers for their comments.

Conflict of interest

The authors declare there is no conflicts of interest.