Approximate optimal iterative approach for a class of oscillatory neuron models

ABSTRACT Control of oscillatory neuron models becomes a growing interest due to the application of implanted stimulus electrodes in mitigating pathological behaviours. We present an approximate optimal iterative control method for minimum-current control of the FitzHugh–Nagumo model with external disturbance. We first focus on revealing optimality conditions based on the Pontryagin's maximum principle and constructing a related nonlinear two-point boundary value problem. Then, we transform the problem into two iterative sequences of linear differential equations. The control law is finally derived and composed of feedback and compensation terms which can be approximated using approximate iterative approach. Simulations illustrate the results.


Introduction
Recently, control of neuron systems by external current stimuli/input has drawn considerable attention due to its applications in deep brain simulation and oscillatory neurocomputers (Dasanayake & Li, 2011;Handa & Sharma, 2016;Osipov, Kurths, and Zhou, 2007;Uhlhaas & Singer, 2006). Many results, concerning the neuron energy efficiency (Niven & Laughlin, 2008;Sengupta, Stemmler, Laughlin, & Niven, 2010), have been developed related to optimal electrical stimulation in single neuron for different control performances, such as minimum energy and spike times. In the field of theoretical neuroscience and automatic control, many control theoretic methods have been employed for optimal controls of spiking neurons (Dasanayake & Li, 2015;Lou & Swamy, 2017;Moehlis, Shea-Brown, & Rabitz, 2006;Nabi, Stigen, Moehlis, & Netoff, 2013;Wilson, Holt, Netoff, & Moehlis, 2015). In particular, in the context of phasereduced neuron models, a constrained magnitude input current as the stimulus was applied to cause a neuron to spike with a maximal firing rate in Moehlis et al. (2006). The maximum principle was used to derive chargebalanced minimum-power controls for a general phase model of neuron oscillators (Dasanayake & Li, 2015). In Wilson et al. (2015), the authors developed an energy optimal control strategy to entrain heterogeneous noisy neurons and apply it to numerical models of noisy phase oscillators and to vitro hippocampal neurons. In spite of these advances, it has not received much attention until now with few results (Danzl, Hespanha, & Moehlis, 2009; CONTACT X. Lou louxy@jiangnan.edu.cn; Qian Ye qqianye@126.com Nabi, Mirzadeh, Gibou, & Moehlis, 2013;Yi, Wang, Li, Wei, & Deng, 2016) appearing on control of conductancebased neuron models which are more intricate than phase-reduced models. In particular,  designed an event-based desynchronizing control stimulus for a network of pathologically synchronized coupled neurons and employed a Hamilton-Jacobi-Bellman (HJB) approach to transform the optimal control problem into solving numerical solution of an HJB partial differential equation by an level set method. By using similar HJB framework, Danzl et al. (2009) proposed an eventbased minimum-time optimal control for oscillatory neurons and the results were extended to a network of globally coupled neurons. From the view point of extracellular electric field stimulus, the optimization of stimulus energy for a reduced two-compartment neuron model was addressed in Yi et al. (2016). Particular motivation for the study on optimal control of oscillatory neuron models using approximate iterative-type approach comes from its board application to approximate solution of optimal control problem for nonlinear systems and the rapid development of computer ability.
As an important spiking neuron model, the FitzHugh-Nagumo model was introduced as a simplified model of the cell membrane (FitzHugh, 1961). During recent years, many approaches have been developed for the case of the FitzHugh-Nagumo model with or without the influence of random perturbations, see Barbu, Cordoni, and Di Persio (2016) and Uzunca, Küçükseyhan, Yucel, and Karasözen (2017), where the model is described by partial differential equations with reaction-diffusion terms. In recent years, the study of control in the class of FitzHugh-Nagumo models has attracted some attention, including (Chavarette, Balthazar, Peruzzi, & Rafikov, 2009;Ciszak et al., 2017;Danzl et al., 2009;FitzHugh, 1961;Hilhorst & Rybka, 2010;Liu & Wang, 2017;Sun & Yang, 2016;Uzunca, Küçükseyhan, Yücel, & Karasözen, 2017;Yi et al., 2016;Zheng, 2015Zheng, , 2016. In particular, an eventbased control approach for reduced Hodgkin-Huxley model and FitzHugh-Nagumo model was developed in Danzl et al. (2009) by means of the Hamilton-Jacobi-Bellman approach. The FitzHugh-Nagumo type reactiondiffusion system is also a generic model in physiology, describing complex wave phenomena in excitable or oscillatory media (Uzunca et al., 2017). Optimal control of semi-linear reaction-diffusion equations is an active research field with many applications in controlling pattern formation (Zheng, 2015), boundary control matched disturbance (Liu & Wang, 2017) and physical processes in motion of turbulence (Zheng, 2016), to give a few examples. Although some existing works addressed optimal control or stabilization problems for FitzHugh-Nagumo type reaction-diffusion systems (Hilhorst & Rybka, 2010;Sun & Yang, 2016;Uzunca et al., 2017), in this paper, we shall mainly focus on the FitzHugh-Nagumo model described by a second-order original differential equation as in Chavarette et al. (2009) andCiszak et al. (2017). By means of Sinha's theory and a linear feedback control strategy, Chavarette et al. (2009) reduced the oscillatory movement of a Fitzhugh-Nagumo system to a stable equilibrium point, but the designed optimal controller depended on a limited matrix G. In Ciszak et al. (2017), Ciszak et al. addressed several dynamical regimes of a FitzHugh-Nagumo model driven by a feedback control composed of a high-pass or an all-pass filter, and explored the relationship between complex dynamical regimes and control parameters (such as the external driving amplitude A, time constant τ and coupling strength coefficient α). We believe that efficient optimal control approaches for FitzHugh-Nagumo systems should play a more prominent role in analysis and control of conductance-based neuron models. Existing approximation iteration approaches for optimal control of nonlinear systems, such as Galerkin approximation (Beard, Saridis, & Wen, 1997;Kim, Kim, & Lim, 2005) and Riccati-equation approximation sequence method (Banks, 1986;Çimen & Banks, 2004), have limitations in the lack of convergence or in high computation complexity.
In this paper, inspired from Çimen and Banks (2004), Sun and Fan (2008) and Gao, Zhang, and Zhang (2008), we aim at developing an approximate iterative approach (AIA) that can be applied to optimal control of a dimensionless FitzHugh-Nagumo model with external disturbance so as to guarantee the minimum-current suboptimal performance. The optimal control law is composed of an exact linear feedback and a nonlinear compensation term and the convergence of the developed approximate iterative algorithm is proved. Since the proposed AIA is related to a pair of vector iterations and the solution of the optimal control is transformed into a limit process of the nonlinear compensated vector sequence, the computational complexity can be greatly reduced. By simulation, we compare the results with those obtained by an orthogonal collocation method and a state-dependent Riccati equation approach, which illustrates that the results using the proposed approach are comparable to those obtained by other methods.

FitzHugh-Nagumo model
As is well-known, the Hodgkin-Huxley model is very important in describing the transmission of an action potential through a cell membrane (Reddy & Murthy, 2013). However, due to the large number of variables, the phase space dynamics of the equation is hard to be visualized. The FitzHugh-Nagumo model, as a simplification of the Hodgkin-Huxley model, was introduced in FitzHugh (1961) and experimentally verified by Nagumo, Arimoto, and Yoshizawa (1962) using electrical circuits. Here, we consider the FitzHugh-Nagumo neural model (Danzl et al., 2009;Keener & Sneyd, 2008) where v is the voltage potential of the neuron membrane, ω is the inactivation of the sodium channels, u(t) represents the input current, δ acts as the singular perturbation parameter and v evolves on a much faster time scale than ω. The neural system can exhibit periodic spiking dynamics, while neurons consisting of such coupled models may be firing in synchronism and lead to pathological behaviours. An efficient way to desynchronize is to use the phaseless set of system states which is extremely sensitive to background noise (Danzl et al., 2009). For the FitzHugh-Nagumo model (1), the origin is the phaseless set contained by the basin of attraction of the periodic orbit (Danzl et al., 2009), and thus our objective is to construct a minimum-current optimal control law to reach the origin.

Minimum-current optimal control
In this section, we introduce the minimum-current optimal control problem of a stochastic version of the FitzHugh-Nagumo model. The dynamics of the neuron model (1) with external disturbance can be written aṡ where x = [v, ω] ∈ R 2 is the state vector of the system, u(t) is the external input injected to the neuron, ξ(t) is an additional disturbance (which may be sinusoidal disturbance, periodic disturbance, persistent disturbance (Gao, 2008), uniform noise or white Gaussian noise) with In practice, for treatment of pathological behaviours, such as Parkinson's disease (Rubchinsky, Park, & Worth, 2012) and obstructive sleep apnea, electrical stimuli in a deep brain stimulation setting is an efficient way with a requirement of clinically reducing stimulation energy. Motivated by this, the interest of this work is in developing an approximate optimal iterative approach for the neuron model when controlling neurons to achieve certain desired behaviour or driving the neurons to their phaseless sets. Considering the neuron system (2), the objective is to find the minimum-current input u * (t) which would drive the states of the neuron to the phaseless set in some prespecified length of time [t 0 , t f ]. This is achieved by minimizing an objective function where x f is the desired terminal state, Q, Q f ∈ R 2×2 are positive semi-definite matrices and R > 0 is a positive constant. These parameters reflect the relative values of the weights determining the relative importance of minimizing each term; e.g., a high R relative to Q and Q f would emphasize minimization of the input current.
In practice, if one extremely needs to achieve certain desired behaviour or drive the neurons to their phaseless sets, one could set high weight coefficients Q and Q f . Extending the optimal control concepts to the neuron systems is motivated by control applications where a high-power stimulus is harmful to the biological tissues and low-power electrical stimulus plays an important role in reducing the power consumption in a neurological implant. Consider the Hamiltonian for the system (2) given by where λ is the co-state variable vector.
According to the Pontryagin's maximum principle, the optimality conditions are obtained by the following nonlinear two-point boundary value problem (BVP) , and the optimal control law is given by Generally, it is hard to analytically solve the nonlinear two-point BVP (5) except for a few simple cases. In the next section, we introduce a transformation of the twopoint BVP and an approximation iterative approach to overcome this difficulty.

Approximation optimal iterative design process
By (5), let us construct the set of two-point BVPṡ and the corresponding optimal control sequence where x (0) (t) denotes the initial states during [t 0 , t f ] before the first iteration, which is used to calculate f (x (k−1) ) and G(x (k−1) ) when k = 1.
Note that ξ(t) remains unchanged during each iteration as it is an additional disturbance. The following theorem shows the relationship between the kth optimal solution and the optimal state trajectory x * (t) and optimal control law u * (t). (7) and (8), respectively. Then, {x (k) (t)} ∞ 1 and {u (k) (t)} ∞ 1 uniformly converge to the optimal state trajectory x * (t) and the optimal control law u * (t), respectively, for the system (2) with the quadratic objective function (3).
Proof: When k = 1, it follows from (7) that for all t ∈ [t 0 , t f ] For k = 0, 1, . . ., let us define Then, (9) can be written aṡ Since (9) is a nonlinear two-point BVP, let us denote where w (1) (t) is an adjoint vector introduced to compensate the nonlinearity of system and P (1) (t) is the unique positive semi-definite solution to a Riccati matrix differential equation that will specified later. Then, substituting (12) into (11a) giveṡ Calculating the derivative of both sides of (12) and using (13) yieldṡ On the other hand, substituting (12) into (11b) giveṡ Using the equivalence between (14) and (15), we obtain the Riccati matrix differential equation and the adjoint vector differential equation

−ẇ (k) (t) = −P (k) (t)Sw (k) (t) + P (k) (t)Dξ(t),
Substituting (20) into (8), we obtain the kth optimal control law In addition, substituting (20) into (19a), we can get the kth optimal closed-loop systeṁ Using the Cauchy sequence theorem and transfer matrix method, one can readily obtain that solutions to the sequence of (24) uniformly converge to the solution of the following equation and solutions to the sequence of (26) uniformly converge to the solution of the following equatioṅ In addition, using the fact λ(t) = P(t)x(t) + w(t), since the control sequence {u (k) (t)} ∞ 1 depends on the state sequence {x (k) (t)} ∞ 1 and the adjoint vector sequence {w (k) (t)} ∞ 1 , then, {u (k) (t)} ∞ 1 also uniformly converges to the optimal control law u * (t), that is, 1 also uniformly converges to the optimal state trajectory x * (t) of the system. The proof is complete.
Theorem 4.1 proves the convergence of the state sequence and the control sequence depending on the adjoint vector sequence, and the proof of Theorem 4.1 reveals the design process of AIA. Note that the closedloop compensation term in the optimal control law (25) related to the adjoint vector will increase the robustness of the system as the disturbance effect is reflected in the adjoint vector. However, it is impractical to use the sequences as it contains infinite series. To apply the AIA in practice, we propose the following efficient algorithm in Algorithm 1 for solving the minimum-current optimal control problem. In this algorithm, the Mth iteration suboptimal control law is employed by replacing ∞ with a finite positive integer M in (30), i.e.
where M is determined by a concrete error criterion.
Step 2. Calculate G (0) (t) and f (0) (t) by means of the definitions Step 3. Obtain the positive semidefinite matrix P (k) (t) by solving the Riccati matrix differential equation Step 4. Obtain w (k) (t) by solving the adjoint vector differential equation and calculate x (k) (t) from the above equation.
Step 6. Set M = k and calculate u (M) (t) from the control law Step 7. Compute G (k) (t) and f (k) (t) by means of the defini- Step 8. Calculate J (M) given by Step 9. If the following inequality holds for the constant ε > 0 (which represents the tolerance error bound), go to step 10; else let k = k + 1 and go to Step 3.
Step 10. Stop the algorithm; obtain x (M) (t) and u (M) (t) which are accurate enough.
Remark 4.1: It is shown from the Practical Algorithm, during the iterative procedure, the parameter Q affects the value P (k) (t) in the Riccati matrix differential equation, the parameter Q f determines the boundary values P (k) (t f ) and w (k) (t f ), and the parameter R affects the state sequence x (k) (t) and further influences the optimal control law.

Remark 4.2:
To compare with existing approximation iteration approaches for optimal control of nonlinear systems, we give a brief comparative analysis as follows.
• The idea of Galerkin approximation method (Beard et al., 1997;Kim et al., 2005) is to directly transform the HJB equation into a sequence of nonlinear vector differential equations and achieve suboptimal control by solving the limit of the sequence. Although the iterative process is fast, it takes a long time to calculate the inner product related to a cost-to-go function for iteration. In addition, this method depends on the selection of initial values and the iterative process will be slow and even divergent without appropriate selection of initial values. However, such a method is very suitable for solving optimal control problems of nonlinear systems with zero output. • The power series expansion method (Cebubar & Costanza, 1984) constructs a control law by using a power series and obtains the approximate optimal solution in the form of sequence. However, it requires the system is analytical to the state vector so that the power series expansion can be carried out which is unrealistic in many practical systems. • In the state-dependent Riccati equation (SDRE) approach (Pittner & Simaan, 2006), the Riccati matrix differential equation (RMDE) with state vector related to P(x) is transformed into an iterative form of constant-parameter RMDE and the solution sequence or series of the RMDE is used to approximate the solution of SDRE. This approach avoids solving the nonlinear HJB equation itself, but it depends on the approximate solution of P(x) and needs iterative solution of a series of matrix differential equations. Therefore, the computational requirements by using this approach are large. Such an issue also exits in the Riccati-equation approximation sequence method (Banks, 1986;Çimen & Banks, 2004), although this method has a faster convergence speed when the dimension of the system is not very large. • The AIA in Gao (2008) is only applicable to nonlinear systems that can be decomposed into independent linear state parts and bilinear input parts. It has the advantage of simple calculation, but it is sensitive to the external weak disturbance and the parameter perturbations within the system. • In contrast to above approximate methods, the proposed AIA avoids solving directly the nonlinear HJB equation or SDRE and provides a convenient way to obtain the approximate solutions for the suboptimal control problems of oscillatory neuron models. Since the affine nonlinear term G(x) is treated as the 'parameter perturbation' of a linear system, the proportion of the closed-loop feedback term has a much higher effect than the case of treating the nonlinear term in the system as an additional disturbance. Therefore, it is robust to the parameter perturbations within the system. • However, one limitation of applying the AIA is that it only causes great positive effects for nonlinear systems that can be decomposed into independent linear parts (that is, Ax(t)) and state-dependent linear parts (that is, G(x(t))x(t)). Hence, it is not necessary to apply the AIA in neuron models with highly nonlinear functions, such as the Hodgkin-Huxley model. It remains an open problem in our further work.

Simulations
For purpose of illustration, let us consider the FitzHugh-Nagumo model (2) with the objective being to control the states of the neuron system to the phaseless point x f = [0, 0] while maintaining the performance given by (3), where, without loss of generality, the parameters are chosen as t 0 = 0, t f = 1, Q f = 100I, Q = 0 and R = 2. Similar to the FitzHugh-Nagumo model in Lindner, García-Ojalvo, Neiman, and Schimansky-Geier (2004), we consider the external disturbance Dξ(t) driven by the white Gaussian noise with intensity d = 0.1. For simplicity, take δ = 1 to eliminate input and voltage scaling. Following the procedure given in Practical Algorithm, we can obtain the suboptimal control law for the neuron model (2) with initial values v 0 = 1, ω 0 = −1 and the objective function index (3). The objective function values at different iterations are listed in Table 1. It can be found from Table 1 that, the greater the iteration gets, the higher the control precision will be. However, when the iteration increases up to a high value, the values of the objective function tend to a determinate constant.
To obtain an accurate enough suboptimal trajectory, we set the tolerance error bound ε = 0.005. In this case, convergence is achieved after 4 iterations, i.e. |(J (4) − J (3) )/J (4) | = 0.0035 < 0.005 and a minimum of J (4) = 2.9510 is obtained. It is worth to point out that, in the proposed algorithm, only a few iteration times are required in order to get a suboptimal control law. In order to verify whether the derived control achieves the objectives or not, it is applied to the neuron system to obtain controlled state trajectories. We compare the results of the AIA with the state trajectories and the optimal control trajectories obtained by using the SDRE approach (Pittner & Simaan, 2006) and the Gauss pseudospectral method (GPM) (Rao et al., 2010). As an orthogonal collocation method, GPM is an efficient optimization approach based  on approximating the state and control trajectories using interpolating polynomials and has been used extensively throughout the world both in academia and industry https://en.wikipedia.org/wiki/GPOPS-II (Huntington & Rao, 2008). Figures 1 and 2 provide simulation results for minimum-current optimal control of the FitzHugh-Nagumo model using three different methods. Figure 1 shows the state trajectories of the solutions. It is seen from the terminal values of w that, the SDRE approach and our AIA perform better than the GPM, that is, the phaseless point x f are achieved by the SDRE approach and our AIA. This will be further illustrated by the error metric E in Table 2. Figure 2 shows the control trajectories of the solutions. It can be seen that among the three methods, the GPM needs the least energy. However, the AIA and the SDRE approach are comparable with the GPM for the needed energy. This is quantatively confirmed by the computational results listed in Table 2, where the iteration length, the computational time (CT) and the objective function value J for all methods are given. The error between the real and the desired terminal states, defined as E = (v(t f ) − v f ) 2 + (ω(t f ) − ω f ) 2 , is also given in Table 2. Obviously, the proposed AIA requires the least iterations and least computational time with tolerable control performance and good accuracy on the desired state. The SDRE approach takes the longest computational time although it achieves the best accuracy on the desired state.

Conclusions
In this paper, we have studied the minimum-current optimal control of the FitzHugh-Nagumo model with external disturbance. By utilizing the Pontryagin's maximum principle, the optimal control problem has been transformed into solving a nonlinear two-point BVP. By solving the two-point BVP using iterative sequences of linear differential equations, the suboptimal control law has been designed through an approximation optimal iterative process. We have carried out simulations using the proposed method for the FitzHugh-Nagumo model. Comparisons of the results by using the SDRE approach, the GPM and the proposed AIA have also been revealed. An extension effort to more complicated neural systems such as the Hodgkin-Huxley model perturbed by general Lévy type noise is possible, which is left for future work.