Sliding mode control of the Hodgkin-Huxley mathematical model

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.


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 ε (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 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 N a , g l are the maximum conductances of these ions, V K , V N a , V l are the constant equilibrium potentials for these ions, C M is the membrane capacitance, δ = a 2R2 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) ∈ Q := (0, T ) × (0, L), where x represents the longitudinal distance along the axon and t is time.
In (1.2)-(1.4) α n , α m , α h , β n , β m , β h are nonlinear functions of v defined as indicated, e.g., in [5], [12], [13], namely The values δ, C M , g K , g N a , g l , C M are positive numbers and V K , V N a , 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 * at a finite time T * and to keep this value for t ≥ T * . The other state variables n, m, h will have after T * an evolution governed by their equations in which v takes the value v * . 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 * 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 where the symbol sign denotes the multivalued function and ρ is a positive constant. We rewrite (1.1)-(1.4) in the following form The system is completed by homogeneous Neumann boundary conditions for v, ∂v ∂x (t, 0) = ∂v ∂x (t, L) = 0, for t ∈ (0, T ), (1.11) since the membrane potential does not have a flux across the ends of the fiber, and by initial conditions 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 If z ∈ L ∞ (X) the notation z ∞ will stand for z L ∞ (X) , where X can be Ω, 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 dv dt , ∂v ∂x , ∂ 2 v ∂x 2 and similarly for the other functions.

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 The desired final value to be obtained is the time and space dependent function v * . 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: 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 ∞ (Q) and so the integrals containing these functions make sense.
and consider v 0 ∈ V, (2.12) (2.14) 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 ε be positive and introduce the Yosida approximation of the sign operator, and the approximating system Let R be a positive value, which will be later specified, and let us introduce the set . We shall apply the Banach fixed point theorem in M.
We fix (v, w) ∈ M and consider the system Let us set We note that f iM depend on w M while h iR depend on R, i = 1, 2 and take where C is a constant depending on the problem parameters and the initial datum for v ε and will be given below. Next, we define Ψ : Moreover, by (2.21) we see that In order to deal with the parabolic problem (2.20), (2.22), (2.23) we introduce the linear time dependent operator and write the equivalent Cauchy problem The operator A(t) has the properties and so by the Lions theorem (see [21], p. 162), the Cauchy problem has a unique solution y ε ∈ The solution satisfies a first estimate, obtained by testing (2.29) by y ε (t) in H and then integrating over (0, t) where δ 1 = min{1, 2a, 2δ}.
We calculate a second estimate, by multiplying formally (2.29) in H by −(y ε ) xx (t) and then integrating over (0, t). We get The latter together with (2.30) provides where C 2 is given by T and δ 2 = min{1, δ}. Recalling (2.25) we deduce that Next, by (2.29) we calculate We also recall that, by (2.28), and with homogeneous Neumann boundary conditions for (y ε1 − y ε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, where the constant C ε depends on ε. We multiply (2.36) scalarly in H by (y ε1 − y ε2 ) and (2.37) by (z ε1 − z ε2 ). We sum up the resulting equations and then we integrate over (0, t). We get Taking the supremum for t ∈ [0, T ], and choosing γ large enough such that 2γ > C ε , we obtain with C εB = C ε /(2γ) < 1, so that Ψ turns out to be a contraction and to have a unique fixed point, This implies that the pair (v ε , w ε ) = (y ε , z ε ) 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 ε ) we have By the Lions-Aubin lemma (see e.g., [21], p. 58) we get implying v ε → v and w ε → w a.e. on Q. By the continuity of f i and h i and by the Lebesgue dominated convergence theorem we also get Also it follows that sign ε (v ε − v * ) → ζ weak-star in L ∞ (Q) and since sign is weakly-strongly closed we get ζ ∈ sign (v − v * ) a.e. (t, x) ∈ Q, (see e.g., [1], p. 38, Proposition 2.2). Now, we consider the weak formulation of (2.16) with ζ ε = sign ε (v ε − v * ) and pass to the limit as ε goes to zero, obtaining (2.9). To this end we took into account that because v ε → v and f i (w ε ) → f i (w) strongly in L 2 (0, T ; H) and f 1 (w)ψ ∈ L 2 (0, T ; H).
Passing to the limit in the weak formulation on (2.17) and taking into account that )ψdxdt → 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
Then, by integrating on (0, L) we get Since w ε (t) → w(t) strongly in H uniformly in t we have that w verifies and it is clear that w ∈ 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 , x) ∈ Q) and the equations corresponding to w 1 and w 2 , Let us multiply the first difference by v 1 − v 2 and the second by w 1 − w 2 , integrate over (0, t) × (0, L) and sum the resulting equations. After similar calculations as before, we get H ds which yields, by Gronwall's lemma, that v 1 (t) = v 2 (t) and w 1 (t) = w 2 (t), for all t ∈ [0, T ]. This proves the solution uniqueness and ends the proof.
We prove now the occurence of the sliding mode at a finite time T * .

Theorem 2.2. Let
and let Then, for T * ∈ [0, T ] defined as Proof. We shall compare the solution to (2.1) with the solution to the system where (·) + is the positive part. Moreover, it can be noticed that q(T * ) = 0 where T * is given by (2.42).
Observe that the function q is positive and decreasing, |q(t)| ≤ |q(0)| for t < T * , it reaches the value zero at T * and remains zero after T * . It is clear that due to the choice (2.41) we have T * < T. We denote p = v − v * and consider the system by (2.40), where ζ p ∈ sign p and ζ q ∈ sign q. We took into account that p(0) − q(0) = v 0 − v * − v 0 − v * ∞ ≤ 0 and so (p(0) − q(0)) + = 0. From here it follows that p(t) ≤ q(t) for all t ∈ [0, T ]. Now, we add the equations for p and q and multiply their sum by −(p + z) − and integrate over (0, t) × (0, L). Taking into account that signz = −sign(−z) we obtain implying that p(t) ≥ −q(t) for all t ∈ [0, T ]. Finally we have obtained that |p(t)| = |v(t) − v * (t)| ≤ q(t), so that v(t) − v * (t) = 0 for t ≥ T * , 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.

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 k i , 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 , r, r 1 , r 2 , r 3 ∈ R, |r i | ≤ M, |r i | ≤ M, there exist L fi (M ) and L h k i (M ) positive, such that (ii) 1 there exists a > 0 such that together with the initial conditions (1.12).

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) = ρ sign ε (v(t, x) − v * (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 ε.
We considered the domain Q = [0, T ] × [0, L] with L = 1 and T ∈ {100, 200, 400} and the approximation of the multivalued function sign given by 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]).
We used the following stopping criteria of the algorithm: that is, the number of iterations exceeds N iter(= 100); eT ime > N T ime, that is, the elapsed time exceeds the limit N T ime(= 900CP U ).
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.
-Evaluate, in the grid points meshgrid the functions v 0 , n 0 , m 0 , h 0 and v * .
-Compute v iter using the pdepe solver, with the initial datum v 0 and the boundary data.
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 -eT ime = 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, n 0 = 0.45, m 0 = 0.03, h 0 = 0.397.
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]: The values (4.1) are the values of membrane channel conductance which, for ρ = 0, do not lead to an unstable membrane potential response (see [12]). This situation is illustrated in Fig. 1 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 * (t, x) = 0.5 sin 4 π t + 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 ρ = 20 and it is illustrated in Fig. 5. All the other parameters are those from (4.1). The final Fig. 6 shows the evolution of the system towards the sliding stabilization when v 0 = 0.5 sin(4πx)+0.6 and g K = 36. In this case the graphics v(t, x f ixed ) differ when modifying x f ixed ∈ [0, 1] due to the space dependence of v 0 . In all situations starting from a constant v 0 we observe that the graphics v(t, x f ixed ) with x f ixed ∈ [0, 1] are the same, since the system is invariant to the translation x → 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.