Sliding mode control of the Hodgkin-Huxley mathematical model

TL;DR

This paper develops a sliding mode control strategy for the Hodgkin-Huxley neuronal model, using an external current to steer the membrane potential to a desired state and maintain it there, supported by mathematical proofs and simulations.

Contribution

It introduces a novel control approach for the Hodgkin-Huxley model employing relay-based feedback to achieve finite-time reaching and sliding on a target manifold.

Findings

Proves existence of sliding mode in the Hodgkin-Huxley system.

Provides a method to determine the time to reach the target manifold.

Demonstrates effectiveness through numerical simulations.

Abstract

In this paper we deal with a feedback control design for the action potential of a neuronal membrane in relation with the non-linear dynamics of the Hodgkin-Huxley mathematical model. More exactly, by using an external current as a control expressed by a relay graph in the equation of the potential, we aim at forcing it to reach a certain manifold in finite time and to slide on it after that. From the mathematical point of view we solve a system involving a parabolic differential inclusion and three nonlinear differential equations via an approximating technique and a fixed point result. The existence of the sliding mode and the determination of the time at which the potential reaches the prescribed manifold is proved by a maximum principle argument. Numerical simulations are presented.

Equations

$$C_M\frac{dv}{dt}=\delta \frac{{\partial }^{2}v}{\partial x^2}-g_K{n}^4(v-V_K)-g_{Na}{m}^{3}h(v-V_{Na})-g_l(v-V_l)+I_C,\text{mbox}inQ,$$

$$\frac{dn}{dt}={\alpha }_n(v)(1-n)-{\beta }_n(v)n,\text{mbox}inQ,$$

$$\frac{dm}{dt}={\alpha }_m(v)(1-m)-{\beta }_m(v)m,\text{mbox}inQ,$$

$$\frac{dh}{dt}={\alpha }_h(v)(1-h)-{\beta }_h(v)h,\text{mbox}inQ,$$

$${\alpha }_n(v)$$

$${\alpha }_m(v)$$

$${\alpha }_h(v)$$

$$I_C(t,x)=-\rho \text{mbox}sign(v(t,x)-v^{\ast }(t,x)),$$

$$\mbox{sign\thinspace}r=\left\{\begin{tabular}[]{ll}$1,$&$r>0$\\ $[-1,1],$&$r=0$\\ $-1$,&$r<0$\end{tabular}\right.$$

$$C_M\frac{dv}{dt}=\delta \frac{{\partial }^{2}v}{\partial x^2}-f_1(n,m,h)v+f_2(n,m,h)+I_C,\text{mbox}inQ$$

$$\frac{dn}{dt}$$

$$\frac{dm}{dt}$$

$$\frac{dh}{dt}$$

$$f_1(n,m,h)$$

$$f_2(n,m,h)$$

$$h_1^n(v)$$

$$h_2^n(v)$$

$$\frac{\partial v}{\partial x}(t,0)=\frac{\partial v}{\partial x}(t,L)=0,\text{mbox}fort\in (0,T),$$

$$v(0,x)=v_0,\text{mbox}n(0,x)=n_0,\text{mbox}m(0,x)=m_0,\text{mbox}h(0,x)=h_0,x\in (0,L).$$

$$V=H^1(0,L)\subset H=L^2(0,L)\subset V^{\prime }=(H^1(0,L){)}^{\prime }$$

$$W=\{y\in H^2(0,L);\text{mbox}{y}_x(0)=y_x(L)=0\}.$$

$$v_t-\delta v_{xx}+f_1(w)v+\rho \text{mbox}sign(v-v^{\ast })\ni f_2(w),\text{mbox}inQ,$$

$$w_t=-h_1(v)w+h_2(v),\text{mbox}inQ,$$

$$v_x(t,0)=v_x(t,L)=0,\text{mbox}in(0,T),$$

$$v(0,x)=v_0,\text{mbox}w(0,x)=w_0,\text{mbox}in(0,L).$$

$$\left|f_i(r)-f_i(\overline{r})\right|\le L_{f_i}(M)\left|r-\overline{r}\right|,\text{mbox}\left|h_i(r)-h_i(\overline{r})\right|\le L_{h_i}(M)\left|r-\overline{r}\right|;$$

$$0<a\le f_1(r),\text{mbox}0<h_i(r),\text{mbox}i=1,2,\text{mbox}forallr\in \mathbb{R};$$

$$v_0\in L^2(0,L),\text{mbox}{w}_0\in L^2(0,L);$$

$$v^{\ast }\in W^{1,\infty }(0,T;H)\cap L^{\infty }(0,T;W).$$

$$v\in C([0,T];H)\cap L^2(0,T;V)\cap W^{1,2}(0,T;H)\cap L^{\infty }(Q),\text{mbox}w\in L^{\infty }(Q)\cap W^{1,2}(0,T;H)$$

$${\int }_0^T{\int }_0^L(v_t\psi +\delta \nabla v\cdot \nabla \psi +f_1(w)v\psi +\rho \zeta \psi )dxdt$$

$$\zeta (t,x)$$

$${\int }_0^T{\int }_0^L(w_t\varphi +h_1(v)w\varphi )dxdt={\int }_0^T{\int }_0^L{h}_2(v)\varphi dxdt,\text{mbox}forall\varphi \in L^2(0,T;H),$$

$$\underset{r\in \mathbb{R}}{\sup }\frac{h_2(r)}{h_1(r)}=:w_M\in \mathbb{R},\text{mbox}$$

Full text

Sliding mode control of the Hodgkin–Huxley mathematical model

Cecilia Cavaterra1, Denis Enăchescu2, Gabriela Marinoschi2,3

1Dipartimento di Matematica “F. Enriques”, Università degli Studi di Milano

Via C. Saldini 50, 20133 Milano, Italy

and

Istituto di Matematica Applicata e Tecnologie Informatiche “E. Magenes”, CNR

Via Ferrata 1, 27100 Pavia, Italy

cecilia.cavaterra@unimi.it

2“Gheorghe Mihoc-Caius Iacob” Institute of Mathematical Statistics and Applied Mathematics of the Romanian Academy,

Calea 13 Septembrie 13, Bucharest, Romania

denaches@fmi.unibuc.ro

gabriela.marinoschi@acad.ro

3Research Group of the Project PN-III-P4-ID-PCE-2016-0372,

Simion Stoilow Institute of Mathematics

of the Romanian Academy, Bucharest, Romania

Abstract. In this paper we deal with a feedback control design for the action potential of a neuronal membrane in relation with the non-linear dynamics of the Hodgkin-Huxley mathematical model. More exactly, by using an external current as a control expressed by a relay graph in the equation of the potential, we aim at forcing it to reach a certain manifold in finite time and to slide on it after that. From the mathematical point of view we solve a system involving a parabolic differential inclusion and three nonlinear differential equations via an approximating technique and a fixed point result. The existence of the sliding mode and the determination of the time at which the potential reaches the prescribed manifold is proved by a maximum principle argument. Numerical simulations are presented.

Keywords: Hodgkin–Huxley model, sliding mode control, feedback stabilization, nonlinear parabolic equations, reaction-diffusion systems

MSC2010 Subject classification: 35K55, 35K57, 35Q92, 93B52, 92C30

1 Introduction

The Hodgkin-Huxley (HH) model is the first complete mathematical model of neuronal membrane dynamics explaining the ionic mechanisms determining the initiation and propagation of action potentials in the squid giant axon. It was successfully established in [20] and since then it has become a prototype model for all kinds of excitable cells, such as neurons and cardiac myocytes. Detailed explanations of the biophysical process illustrated by HH model can be found, e.g., in [3], [5], besides the original work [20]. In [5] an analysis of the non-linear dynamics in the Hodgkin-Huxley mathematical model showing the existence of transient chaotic solutions in the model with their original parameters, combined with the presentation of some modifications in the dynamic system in order to become more realistic, has been done.

Many papers have been devoted to the mathematical analysis of this system which exhibits a very complicated behavior. We confine ourselves to mention some fundamental mathematical works on the traditional Hodgkin-Huxley equations: [10], [14]-[17], [22], [18]. In the last one the existence of a unique classical solution of the Hodgkin-Huxley system was proved. In the paper [4], the authors consider a singular perturbation of the Hodgkin-Huxley system and study the associated dynamical system on a suitable bounded phase space, when the perturbation parameter $\varepsilon$ (i.e., the axon specific inductance) is sufficiently small, proving the existence of bounded absorbing sets, of smooth attracting sets, as well as the existence of a smooth global attractor.

From the mathematical point of view various properties of the dynamics of the Hodgkin-Huxley vector field have been studied. Many studies in the literature reveal bifurcations generated in the HH model, such as Hopf bifurcation, period-double bifurcation and double cycle bifurcation (see e.g., [6] and the references there indicated). The HH model can even exhibit a chaotic regime through a series of bifurcations. The qualitative change of neuronal membrane potential from resting to repetitive spiking, which is a characteristic behavior of this model, is of a particular interest, because abnormal repetitive spiking are proper to several neurological diseases. Consequently, much attention was directed to provide mathematical results aiming to avoid instability around bifurcations or to obtain desired dynamical behaviors which might be of help in the development of the therapies of the diseases. For example, various dynamic feedback control methods have been proposed to control the onset of Hopf bifurcation in HH model, see, e.g., [7] and the references there indicated. We also cite the work [12], where the aim was to develop a novel current control law with the purpose to stop the repetitive firing caused by channel conductance deviations and the work [13], focusing on the simulation of the feedback controlled nerve fiber stimulation where the behavior of the nerve fiber is manipulated by an electrical field generator.

The Hodgkin-Huxley model introduced in [20], p. 522, eq. (29) reads

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $v$ is the electrical potential in the nerve, $n,$ $m$, $h$ are the proportions of the activating molecules of the potassium $(n),$ sodium channels ($m)$ and of the inactivating molecules of the sodium channels $(h)$, respectively, $g_{K},$ $g_{Na},$ $g_{l}$ are the maximum conductances of these ions, $V_{K},$ $V_{Na},$ $V_{l}$ are the constant equilibrium potentials for these ions, $C_{M}$ is the membrane capacitance, $\delta=\frac{a}{2R_{2}}$ is a constant (depending on the fiber radius $a$ and the specific resistance of the axoplasm $R_{2})$ and $I_{C}$ is the applied current. Here $(t,x)\in Q:=(0,T)\times(0,L)$, where $x$ represents the longitudinal distance along the axon and $t$ is time.

In (1.2)-(1.4) $\alpha_{n},$ $\alpha_{m},$ $\alpha_{h},$ $\beta_{n},$ $\beta_{m},$ $\beta_{h}$ are nonlinear functions of $v$ defined as indicated, e.g., in [5], [12], [13], namely

[TABLE]

The values $\delta,$ $C_{M},$ $g_{K},$ $g_{Na},$ $g_{l},$ $C_{M}$ are positive numbers and $V_{K},$ $V_{Na},$ $V_{l}$ are real numbers.

This paper involves a new control approach, the sliding mode control, in order to stabilize the membrane potential to a desired value. Sliding mode control is an efficient tool for the stabilization of continuous or discrete time systems. It consists in finding an appropriate control able to constrain the evolution of the system in such a way to force it to reach a manifold of a lower dimension, called the sliding manifold, in finite time, and to keep it further sliding on this surface. Thus, our purpose is to control the potential $v$ by means of a certain control $I_{C}$ in order to force the potential to reach a prescribed value $v^{\ast}$ at a finite time $T^{\ast}$ and to keep this value for $t\geq T^{\ast}.$ The other state variables $n,m,h$ will have after $T^{\ast}$ an evolution governed by their equations in which $v$ takes the value $v^{\ast}.$ The principal advantage of a sliding mode technique is that after some time the system evolves on a manifold of lower dimension. For recent results regarding sliding mode control for systems of parabolic equations we refer the reader to the papers [2], [8], [9].

The objective is that $v$ reaches a constant value, in particular zero, and to prove in this way the possibility to control the repetitive firing in nerve fibers modeled by the Hodgkin-Huxley system. Even if a constant target might be of main interest, the proof will be developed for a more general case with $v^{\ast}$ dependending on time and space, which allows the target to vary in time, being, for instance, periodic. To this end we propose a relay feedback control of the form

[TABLE]

where the symbol sign denotes the multivalued function

[TABLE]

and $\rho$ is a positive constant.

We rewrite (1.1)-(1.4) in the following form

[TABLE]

[TABLE]

where

[TABLE]

The system is completed by homogeneous Neumann boundary conditions for $v,$

[TABLE]

since the membrane potential does not have a flux across the ends of the fiber, and by initial conditions

[TABLE]

We shall approach this problem in two steps. First, as all equations for the three components $n,$ $m$ and $h$ are similar, we shall consider a reduced system formed only of two equations, one for the potential and the other for only one ionic component, denoted generically by $w.$ This simplification also occurs in the papers Fitzgibbon et al. (see [18], [19]). In Section 2, we shall treat the simplified problem via an approximating method, using a fixed point technique for proving the existence of a solution to the system formed by the equation for the membrane potential, with (1.6) replaced by involving the Yosida approximation and one equation of the form (1.9). Suitable estimates and compactness properties will allow to pass to the limit and to prove an existence result for the non-approximated system in Theorem 2.1. Then, the existence of the sliding mode will be provided in Theorem 2.2 by a comparison argument. In Section 3, we shall extend the result to the complete system (1.8)-(1.9), by observing that it follows as a consequence of the previous results for the simplified system. The paper is concluded by numerical simulations intended to put into evidence the sliding mode behavior of the solution.

Notation. We denote

[TABLE]

where $V\subset H\subset V^{\prime}$ with compact injections. Moreover, we define

[TABLE]

If $z\in L^{\infty}(X)$ the notation $\left\|z\right\|_{\infty}$ will stand for $\left\|z\right\|_{L^{\infty}(X)}$, where $X$ can be $\Omega,$ or $Q.$ We denote by $C,$ $C_{i},$ $i=1,2,...$ some constants depending on problem parameters, sometimes explicitly indicated in the argument. For the sake of simplicity we shall write $v_{t},$ $v_{x},$ $v_{xx}$ instead of $\frac{dv}{dt},$ $\frac{\partial v}{\partial x},$ $\frac{\partial^{2}v}{\partial x^{2}}$ and similarly for the other functions.

2 The simplified system

Let us consider the system for the potential $v$ and the concentration $w,$ coupled with a set of homogeneous Neumann boundary conditions for the potential and of initial data

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The desired final value to be obtained is the time and space dependent function $v^{\ast}$. Here the value $C_{M}$ is considered for simplicity equal to 1.

Taking into account the general considerations presented in the introduction on the expressions of the functions occurring in the Hodgkin-Huxley model, we shall assume the following properties:

• (i)

the functions $f_{i}$ and $h_{i}$, for $i=1,2,$ are locally Lipschitz continuous, that is, for any $M$ positive, and for any $r$, $\overline{r}\in\mathbb{R}$, $\left|r\right|\leq M,$ $\left|\overline{r}\right|\leq M,$ there exist $L_{f_{i}}(M)$ and $L_{h_{i}}(M)$ positive, such that

[TABLE]

• (ii)

there exists $a>0$ such that

[TABLE]

• (iii)

[TABLE]

• (iv)

[TABLE]

Definition 2.1. We call a solution to system (2.1)-(2.4) a pair $(v,w)$

[TABLE]

which satisfies

[TABLE]

[TABLE]

together with the initial conditions (2.4).

We observe that by hypotheses (i) and (2.8) it follows that $h_{i}(v)$ and $f_{1}(w)$ belong to $L^{\infty}(Q)$ and so the integrals containing these functions make sense.

Theorem 2.1. Let (i)-(iv)* hold. Assume that*

[TABLE]

and consider

[TABLE]

[TABLE]

[TABLE]

Then, problem (2.1)-(2.4) has a unique solution, with the further regularity

[TABLE]

**Proof. **We shall consider a regularized problem and prove that it has a unique solution by applying the Banach fixed point theorem. Then, we shall pass to the limit to recover the solution to (2.1)-(2.4).

Let $\varepsilon$ be positive and introduce the Yosida approximation of the sign operator,

[TABLE]

and the approximating system

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Let $R$ be a positive value, which will be later specified, and let us introduce the set

[TABLE]

which obviously is a closed subset of $C([0,T];H)\times C([0,T];H)$. Also, ${\mathcal{M}}$ is a metric space with the metric $d_{\mathcal{M}}((u,w),(\overline{u},\overline{w}))=\left\|u-\overline{u}\right\|_{C([0,T];H)}+\left\|w-\overline{w}\right\|_{C([0,T];H)}.$ We shall apply the Banach fixed point theorem in $\mathcal{M}.$

We fix $(\overline{v},\overline{w})\in\mathcal{M}$ and consider the system

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Since

[TABLE]

with $\overline{v}\in L^{\infty}(0,T;C[0,L])$ (indeed $V\subset C[0,L]$ in a one-dimensional space), we get $f_{i}(\overline{w})\in L^{\infty}(Q)$ and $h_{i}(\overline{v})\in L^{\infty}(Q).$ Let us set

[TABLE]

We note that $f_{iM}$ depend on $w_{M}$ while $h_{iR}$ depend on $R,$ $i=1,2$ and take

[TABLE]

where $\overline{C}$ is a constant depending on the problem parameters and the initial datum for $v_{\varepsilon}$ and will be given below.

Next, we define $\Psi:\mathcal{M}\rightarrow L^{2}(0,T;H)\times L^{2}(0,T;H),$ by $\Psi(\overline{v},\overline{w})=(y_{\varepsilon},z_{\varepsilon})$ the solution to (2.20)-(2.23) and prove further that $\Psi(\mathcal{M})\subset\mathcal{M}$ and that $\Psi$ is a contraction.

By (2.21) we have

[TABLE]

It is immediately seen that $z_{\varepsilon}\in C([0,T];C[0,L])$. Indeed $(t,x)\rightarrow\int_{s}^{t}h_{1}(\overline{v}(\sigma,x))d\sigma$ and $(t,x)\rightarrow\int_{0}^{t}e^{-\int_{s}^{t}h_{1}(\overline{v}(\sigma,x))d\sigma}h_{2}(\overline{v}(\sigma,x))d\sigma$, for $0\leq s<t\leq T$, are continuous on $[0,T]$ and $w_{0}\in C([0,L])$. By (2.6), (2.13), (2.11) it follows that $z_{\varepsilon}(t,x)\geq 0$ and

[TABLE]

hence

[TABLE]

Moreover, by (2.21) we see that

[TABLE]

In order to deal with the parabolic problem (2.20), (2.22), (2.23) we introduce the linear time dependent operator $A(t):V\rightarrow V^{\prime},$

[TABLE]

and write the equivalent Cauchy problem

[TABLE]

The operator $A(t)$ has the properties

[TABLE]

and so by the Lions theorem (see [21], p. 162), the Cauchy problem has a unique solution $y_{\varepsilon}\in W^{1,2}(0,T;V^{\prime})\cap L^{2}(0,T;V)\cap C([0,T];H).$ The solution satisfies a first estimate, obtained by testing (2.29) by $y_{\varepsilon}(t)$ in $H$ and then integrating over $(0,t)$

[TABLE]

where $\delta_{1}=\min\{1,2a,2\delta\}.$

We calculate a second estimate, by multiplying formally (2.29) in $H$ by $-(y_{\varepsilon})_{xx}(t)$ and then integrating over $(0,t).$ We get

[TABLE]

whence

[TABLE]

The latter together with (2.30) provides

[TABLE]

where $\overline{C}^{2}$ is given by

[TABLE]

and $\delta_{2}=\min\{1,\delta\}$. Recalling (2.25) we deduce that

[TABLE]

Next, by (2.29) we calculate

[TABLE]

We also recall that, by (2.28), $w_{\varepsilon}\in C^{1}([0,T];C[0,L])$ and

[TABLE]

Estimates (2.33)-(2.35) and (2.27) ensure that the solution $(y_{\varepsilon},z_{\varepsilon})$ to (2.20)-(2.23) belongs to $\mathcal{M}.$

Now, let us consider two pairs $(v^{1},w^{1}),$ $(v^{2},w^{2})\in\mathcal{M}$, with the same initial data. We denote by $(y_{\varepsilon 1},z_{\varepsilon 1})$ and $(y_{\varepsilon 2},z_{\varepsilon 2})$ the corresponding solutions to (2.20)-(2.23) and we calculate the difference of equations (2.20) and (2.21). Namely, we write

[TABLE]

[TABLE]

with homogeneous Neumann boundary conditions for $(y_{\varepsilon 1}-y_{\varepsilon 2})$ and zero initial data. Relying on the local Lipschitz continuity of $f_{i}$ and $h_{i}$, we perform a few calculations in the right-hand sides of the above equations, denoted $RHS$ and $RHS_{1},$ respectively,

[TABLE]

and

[TABLE]

where the constant $C_{\varepsilon}$ depends on $\varepsilon$. We multiply (2) scalarly in $H$ by $(y_{\varepsilon 1}-y_{\varepsilon 2})$ and (2.37) by $(z_{\varepsilon 1}-z_{\varepsilon 2})$. We sum up the resulting equations and then we integrate over $(0,t).$ We get

[TABLE]

so that, defining

[TABLE]

by Gronwall’s lemma we obtain, for $t\in[0,T],$

[TABLE]

where ${\widetilde{C}_{\varepsilon}}=C(R,w_{M},\rho,h_{1},h_{2},f_{1},f_{2},\varepsilon)$. In order to show that $\Psi$ is a contraction, we introduce further the norm $\left\|q\right\|_{B}=\sup\limits_{t\in[0,T]}\left(e^{-\gamma t}\left\|q(t)\right\|_{H}\right)$ which is equivalent to the standard norm in $C([0,T];H)$. So that, we multiply (2.38) by $e^{-2\gamma t}$ getting

[TABLE]

Taking the supremum for $t\in[0,T],$ and choosing $\gamma$ large enough such that $2\gamma>{\widetilde{C}_{\varepsilon}},$ we obtain

[TABLE]

[TABLE]

with ${\widetilde{C}_{\varepsilon B}}={\widetilde{C}_{\varepsilon}}/(2\gamma)<1,$ so that $\Psi$ turns out to be a contraction and to have a unique fixed point, $\Psi(\overline{v},\overline{w})=(\overline{v},\overline{w})=(y_{\varepsilon},z_{\varepsilon}).$

This implies that the pair $(v_{\varepsilon},w_{\varepsilon})=(y_{\varepsilon},z_{\varepsilon})$ is the unique solution to system (2.16)-(2.19) satisfying the estimates (2.27)-(2.28) and (2.32)-(2.35).

Therefore, along a subsequence (denoted still by ${}_{\varepsilon})$ we have

[TABLE]

By the Lions-Aubin lemma (see e.g., [21], p. 58) we get

[TABLE]

implying $v_{\varepsilon}\rightarrow v$ and $w_{\varepsilon}\rightarrow w$ a.e. on $Q.$ By the continuity of $f_{i}$ and $h_{i}$ and by the Lebesgue dominated convergence theorem we also get

[TABLE]

Also,

[TABLE]

By Arzelà-Ascoli theorem we still obtain

[TABLE]

Also it follows that sign${}_{\varepsilon}(v_{\varepsilon}-v^{\ast})\rightarrow\zeta$ weak-star in $L^{\infty}(Q)$ and since sign is weakly-strongly closed we get

[TABLE]

(see e.g., [1], p. 38, Proposition 2.2).

Now, we consider the weak formulation of (2.16)

[TABLE]

with $\zeta_{\varepsilon}=$ sign${}_{\varepsilon}(v_{\varepsilon}-v^{\ast})$ and pass to the limit as $\varepsilon$ goes to zero, obtaining (2). To this end we took into account that

[TABLE]

because $v_{\varepsilon}\rightarrow v$ and $f_{i}(w_{\varepsilon})\rightarrow f_{i}(w)$ strongly in $L^{2}(0,T;H)$ and $f_{1}(w)\psi\in L^{2}(0,T;H).$

Passing to the limit in the weak formulation on (2.17)

[TABLE]

and taking into account that $\int_{0}^{T}\int_{0}^{L}(h_{1}(v_{\varepsilon})w_{\varepsilon}-h_{1}(v)w))\psi dxdt\rightarrow 0$, in a similar way as above, we get (2.10). These two last equations prove that $(v,w)$ is a solution to (2.1)-(2.4).

Moreover, by straightforward calculations using (2.17) and (2.2) we obtain

[TABLE]

Then, by integrating on $(0,L)$ we get

[TABLE]

Since $w_{\varepsilon}(t)\rightarrow w(t)$ strongly in $H$ uniformly in $t$ we have that $w$ verifies

[TABLE]

and it is clear that $w\in C([0,T];C[0,L]),$ since each term is continuous.

For the uniqueness, let $(v_{1},w_{1}),$ $(v_{2},w_{2})$ be two solutions to (2.1)-(2.4) corresponding to the same initial data. We subtract the equations corresponding to $v_{1}$ and $v_{2},$

[TABLE]

(where $\zeta_{1}\in$ sign $(v_{1}-v^{\ast})$ and $\zeta_{2}\in$ sign $(v_{2}-v^{\ast})$ a.e. $(t,x)\in Q)$ and the equations corresponding to $w_{1}$ and $w_{2},$

[TABLE]

Let us multiply the first difference by $v_{1}-v_{2}$ and the second by $w_{1}-w_{2}$, integrate over $(0,t)\times(0,L)$ and sum the resulting equations. After similar calculations as before, we get

[TABLE]

which yields, by Gronwall’s lemma, that $v_{1}(t)=v_{2}(t)$ and $w_{1}(t)=w_{2}(t),$ for all $t\in[0,T].$ This proves the solution uniqueness and ends the proof. $\square$

We prove now the occurence of the sliding mode at a finite time $T^{\ast}$.

Theorem 2.2. Let

[TABLE]

and let

[TABLE]

Then, for $T^{\ast}\in[0,T]$ defined as

[TABLE]

it holds

[TABLE]

Proof. We shall compare the solution to (2.1) with the solution to the system

[TABLE]

Since $\frac{A}{\rho}<1$ by (2.41), $\frac{A}{\rho}\in$ sign [math] and one can verify that the solution to (2.44) is

---

[TABLE]

where $(\cdot)^{+}$ is the positive part. Moreover, it can be noticed that $q(T^{\ast})=0$ where $T^{\ast}$ is given by (2.42). Observe that the function $q$ is positive and decreasing, $\left|q(t)\right|\leq\left|q(0)\right|$ for $t<T^{\ast},$ it reaches the value zero at $T^{\ast}$ and remains zero after $T^{\ast}.$ It is clear that due to the choice (2.41) we have $T^{\ast}<T.$

We denote $p=v-v^{\ast}$ and consider the system

[TABLE]

[TABLE]

[TABLE]

with homogeneous Neumann boundary conditions both for $p$ and $q.$ Observe that since $q$ depends only on time (cf. (2.46)), then (2.44)-(2.45) is equivalent to (2.48)-(2.49-ii).

We subtract (2.48) from (2.47) and multiply the difference equation scalarly in $H$ by the positive part $(p-q)^{+}$ and integrate over $(0,t).$ By few calculations and majorating the right-hand side of the difference equation, we obtain

[TABLE]

by (2.40), where $\zeta_{p}\in$ sign$\,p$ and $\zeta_{q}\in$ sign$\,q.$ We took into account that $p(0)-q(0)=v_{0}-v^{\ast}-\left\|v_{0}-v^{\ast}\right\|_{\infty}\leq 0$ and so $(p(0)-q(0))^{+}=0.$ From here it follows that $p(t)\leq q(t)$ for all $t\in[0,T].$

Now, we add the equations for $p$ and $q$ and multiply their sum by $-(p+z)^{-}$ and integrate over $(0,t)\times(0,L).$ Taking into account that sign$z=-$sign$(-z)$ we obtain

[TABLE]

where $\zeta_{q}\in$ sign$\,(-q)$. Observe that ($p(0)+q(0))=v_{0}-v^{\ast}+\left\|v_{0}-v^{\ast}\right\|_{\infty}\geq 0$ and so $(p(0)+q(0))^{-}=0.$ Thus,

[TABLE]

implying that $p(t)\geq-q(t)$ for all $t\in[0,T].$

Finally we have obtained that $\left|p(t)\right|=\left|v(t)-v^{\ast}(t)\right|\leq q(t),$ so that $v(t)-v^{\ast}(t)=0$ for $t\geq T^{\ast},$ which yields (2.43), as claimed. Finally, we observe that $A$ is not a sharp value (it could be smaller) but here the objective was to prove its existence. $\square$

3 The complete system

Relying on the results previously obtained we can pass to the complete Hodgkin-Huxley system (1.8)-(1.12) and assume:

• (i)1

the functions $f_{i}$ and $h_{i}^{k},$ $i=1,2,$ $k=n,$ $m,$ $h,$ are locally Lipschitz continuous, that is, for any $M$ positive, and for any $\,r,$ $r_{1}$, $r_{2},$ $r_{3},$ $\overline{r},$ $\overline{r}_{1},$ $\overline{r}_{2},$ $\overline{r}_{3}\in\mathbb{R}$, $\left|r_{i}\right|\leq M,$ $\left|\overline{r}_{i}\right|\leq M,$ there exist $L_{f_{i}}(M)$ and $L_{h_{i}^{k}}(M)$ positive, such that

[TABLE]

• (ii)1

there exists $a>0$ such that

[TABLE]

• (iii)1

[TABLE]

[TABLE]

[TABLE]

• (iv)1

[TABLE]

Definition 3.1. We call a solution to system (1.8)-(1.12) a vector $(v,n,m,h)$

[TABLE]

which satisfies

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

together with the initial conditions (1.12).

Theorem 3.1. Let (1)1-(iv)1* hold. Then, problem* (1.8)-(1.12) has a unique solution, which has the supplementary regularity

[TABLE]

Moreover, if

[TABLE]

with $A$ as in (2.40), *then for * $T^{\ast}\in(0,T)$ defined as

[TABLE]

it holds

[TABLE]

**Proof. In **(1.8)-(1.12) we set

[TABLE]

and so system (1.8)-(1.12) can be written in the form (2.1)-(2.4). Also, we observe that if \overline{w}=\left(\begin{array}[]{c}\overline{n}\\ \overline{m}\\ \overline{h}\end{array}\right) and $\left|n\right|\leq M,$ $\left|\overline{n}\right|\leq M,$ $\left|m\right|\leq M,$ $\left|\overline{m}\right|\leq M,$ $\left|h\right|\leq M,$ $\left|\overline{h}\right|\leq M$, then

[TABLE]

Moreover, $h_{2}(v)$ and each column vector in $h_{1}(v)$ are locally Lipschitz due to the same properties of $h_{i}^{k}(v)$ by (3.1). Then, one can apply Theorem 2.1 and take as set $\mathcal{M}$ the following

[TABLE]

Here, we set $f_{iM}:=f_{i}(0,0,0)+L_{f_{i}}\left(n_{M}+m_{M}+h_{M}\right).$ Then, Theorem 2.2 can be applied to get the result. $\square$

4 Numerical simulations

We present some numerical simulations intended to show the feature of the HH system evolution controlled by the relay controller and to put into evidence the sliding mode behavior.

The numerical simulations have been done for the complete system (1.8)-(1.12) with $I_{C}(t,x)=\rho\,$sign${}_{\varepsilon}(v(t,x)-v^{\ast}(t,x)),$ which was solved by an interactive technique. Thus, the numerical solution is computed for the approximating system, but for simplicity, we shall refer later to $v,$ $n,$ $m,$ $h$ without the subscript $\varepsilon.$

We considered the domain $Q=[0,T]\times[0,L]$ with $L=1$ and $T\in\{100,200,400\}$ and the approximation of the multivalued function sign given by

[TABLE]

with $\varepsilon=10^{-4}.$

To solve the ordinary differential equations, ODE, in $n$, $m$ and $h$ with the corresponding initial conditions we used the ode45 solver from Matlab. The solver is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair ([11]). It is a one-step solver for computing a value $y(t_{n})$ and it needs only the solution at the immediately preceding time point $y(t_{n-1}).$

The numerical solution of the initial-homogeneous Neumann boundary value problem for the 1-D parabolic equation in $v$ was computed with the* pdepe* Matlab solver. The solver discretizes the space using a given $xmesh=[x_{1}<...<x_{maxX}]$ and integrates the resulting ODE to obtain approximate solutions at times specified by a vector of points $tmesh=[t_{1}<...<t_{maxT}]$ for all points in $xmesh$ (see [24]). The time integration is done with ode15s, a variable order Matlab solver based on the numerical differentiation formulas (see [23]).

To discretize the space interval $[0,L]$ we considered $maxX$ steps, with $maxX\in\{5,25\}$, and to discretize the time interval $[0,T]$ we took $maxT$ steps, with $maxT\in\{20,40,200\}$.

We used the following stopping criteria of the algorithm:

[TABLE]

where $v^{iter}$ is the value of $v$ at iteration $iter$;

[TABLE]

that is, the number of iterations exceeds $Niter(=100)$;

[TABLE]

that is, the elapsed time exceeds the limit $NTime(=900CPU)$.

The idea of the algorithm is to solve iteratively, on blocks, the system $(v,n,m,h)$ untill one of the stopping conditions is fulfilled. In our case we consider two blocks, the first containing the PDE in $v$, the second block containing the ODE system $(n,m,h)$.

A solving iteration consists in: computing a numerical solution of $v$ using the previous values of $n$, $m$, $h$ and after, to plug-in the obtained value of $v$ in the right terms of ODE for $n$, $m$, $h$ and solve the equations from the second block.

The pseudo-code of algorithm is:

Step 1. Discretize the domain $Q$, compute the initial and boundary conditions and initialize the solving loop.

• Generate the vectors $tmesh$, $xmesh$ and the matrix $meshgrid=\{(t_{i},x_{j})\}_{i=1,.,maxT}^{j=1,.,maxX}$ using the given $maxT$ and $maxX$.

• Evaluate, in the grid points $meshgrid$ the functions $v_{0},$ $n_{0},$ $m_{0},$ $h_{0}$ and $v^{\ast}$.

• $iter=0.$

• $n^{iter}=n_{0},$ $m^{iter}=m_{0},$ $h^{iter}=h_{0},$ $v^{iter}=v_{0},$ $eTime=0.$

Step 2.

• Compute $f_{1}(n^{iter},m^{iter},h^{iter})$ and $f_{1}(n^{iter},m^{iter},h^{iter}).$

Step 3. Iteration loop.

• $iter=iter+1$.

• Compute $v^{iter}$ using the pdepe solver, with the initial datum $v_{0}$ and the boundary data.

Step 4.

• Evaluate $h_{n}^{1}(v^{iter}),$ $h_{m}^{1}(v^{iter}),h_{h}^{1}(v^{iter})$ and $h_{n}^{2}(v^{iter}),$ $h_{m}^{2}(v^{iter}),h_{h}^{2}(v^{iter}).$

Step 5. For each point in $xmesh$, compute:

• $n^{iter}$ using the ode45 solver, with the initial datum $n_{0}$, for all points in $tmesh$

• $m^{iter}$ using the ode45 solver, with the initial datum $m_{0}$, for all points in $tmesh$

• $h^{iter}$ using the ode45 solver, with the initial datum $h_{0}$, for all points in $tmesh.$

Step 6. Check the stopping criteria.

If one of the three conditions is met the algorithm is finished, otherwise

• $eTime$ = elapsed time from the beginning of the iteration loop.

• Go to Step 2.

The algorithm converges, i.e., the first stopping condition is met, in most cases in a short time (dozen of seconds). The second and third stopping criterion is used in the non-stabilization cases.

The initial condition was selected by assuming that the sodium channel inactivation ratio is higher than that of the activation,

[TABLE]

For the initial $v_{0}$ various values were considered and they are indicated in the figures.

Most part of the parameters used in the computations is the same as in [12]:

[TABLE]

In some figures the values of $\rho,$ $g_{K}$ and $T$ differ from those before and they are specified in the captions.

The graphics plotted in all figures show (from left to write) the time evolution of $v$ at specified fixed $x,$ the surface $v(t,x)$ and the time evolution of the proportions of the activating molecules of the potassium $n$ $(K)$ and natrium $m$ $(Na)$ channels and the proportion of the inactivating molecules of molecules of sodium $h$ $(Iso)$ at specified fixed $x.$

[TABLE]

Fig. 1. Graphics $v(t,0)$ (left), $v(t,x)$ (center), $n$, $m$, $h$ (right) for $v_{0}=4.82$, $v^{*}=0$, $\rho=0$

The values (4.1) are the values of membrane channel conductance which, for $\rho=0,$ do not lead to an unstable membrane potential response (see [12]). This situation is illustrated in Fig. 1. However, for the same values, the computations for $\rho=20$ show in Fig. 2 a quicker stabilization.

[TABLE]

Fig. 2. Graphics $v(t,0)$ (left), $v(t,x)$ (center), $n$, $m$, $h$ (right) for $v_{0}=4.82$, $v^{*}=0$, $\rho=20$

In order to illustrate the theory that allows the solution to reach a periodic sliding mode we present Fig. 3 which describes such a situation for $v^{\ast}(t,x)=0.5\sin\left(\frac{4}{\pi}t\right)+0.6$ and the same values (4.1).

For some deviation in conductance parameters, the situation can change with respect to the first case (see [12]). Thus, for the potassium channel conductance $g_{K}=3.8229,$ the stabilization does not occur and as a matter of fact one observes in Fig. 4 a late firing behavior. The desired stabilization is obtained for a suitable $\rho=20$ and it is illustrated in Fig. 5. All the other parameters are those from (4.1).

[TABLE]

Fig. 3. Graphics $v(t,0)$ (left), $v(t,x)$ (center), $n$, $m$, $h$ (right) for $v_{0}=4.82$, $v^{*}=0.5*sin(4/\pi*t)+0.6$, $\rho=20$

[TABLE]

Fig. 4. Graphics $v(t,0)$ (left), $v(t,x)$ (center),$n$, $m$, $h$ (right) for $v_{0}=4.82$, $v^{*}=0$, $\rho=0$, $g_{K}=3.8229$

[TABLE]

Fig. 5. Graphics $v(t,0)$ (left), $v(t,x)$ (center), $n$, $m$, $h$ (right) for $v_{0}=4.82$, $v^{*}=0$, $\rho=20$, $g_{K}=3.8229$

The final Fig. 6 shows the evolution of the system towards the sliding stabilization when $v_{0}=0.5\sin(4\pi x)+0.6$ and $g_{K}=36.$ In this case the graphics $v(t,x_{fixed})$ differ when modifying $x_{fixed}\in[0,1]$ due to the space dependence of $v_{0}$.

[TABLE]

Fig. 6. Graphics $v(t,0)$ (left), $v(t,x)$ (center), $n$, $m$, $h$ at $x=0.5$ (right) for $v_{0}=0.5sin(4\pi x)+0.6$, $v^{*}=0$, $\delta=50$, $g_{K}=36$, $\rho=50$

In all situations starting from a constant $v_{0}$ we observe that the graphics $v(t,x_{fixed})$ with $x_{fixed}\in[0,1]$ are the same, since the system is invariant to the translation $x\rightarrow x+l$. A different situation can be observed in Fig. 6 when the initial $v_{0}$ depends on $x.$

While the left and center plots in each figure show the evolution of the membrane potential, the right ones put into evidence the play between the other components of the system. The equilibrium potential is determined by gradients of ionic concentration, through the membrane permeability, and also, by the effect of the sodium-potassium transport. There is a concentration of potassium ions inside the cell and a higher concentration of sodium chloride ions in the external part. At their turn, the permeabilities of the membrane to sodium and potassium depend on the membrane potential. The figures on the right show a fast initial inflow of sodium ions and a subsequent outflow of potassium ions, which define the action potential generation that follows the stimulation of the depolarization. The chloride ions do not play their role very well, but they first exhibit an increase.

We proposed a sliding mode control strategy for the Hodgkin–Huxley model, by controlling the equation for the membrane potential by a relay type controller. This permits to reduce the oscillatory movement of the nonlinear Hodgkin–Huxley system to a stable equilibrium point.

Acknowledgments. This research activity has been performed in the framework of the Italian-Romanian project “Control and stabilization problems for phase field and biological systems” of the Italian CNR and the Romanian Academy, 2017-2019 and was partially supported by a grant of Ministry of Research and Innovation, CNCS UEFISCDI, Project Number PN-III-P4-ID-PCE-2016-0372, within PNCDI III, for G.Marinoschi. C. Cavaterra is a member of the Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). D. Enăchescu is a member of the PhD School of Computer Science, University of Bucharest.