Non-linear Model-based Control of Neural Cell Dynamics


 This paper aims to discuss the control of a neuron’s firing. It is desired to control the neuron by injecting a current, whichallows the voltage of the neuron membrane to reach what is called the threshold voltage. When a neuron reaches this voltage,the potassium and sodium ion gates open allowing for the neuron to fire an action potential, observed as a voltage spike.Controlling this effect is helpful for those with certain diseases or disabilities. Four types of controllers are designed andsimulated on both the nonlinear system of the neuron and its linearized form, and all are found to meet specifications.


Introduction
Neurons, the nerve cells in our brains, can have more than millions of connections to other neurons. These neurons are made up of the cell membrane, axons, dendrites, and synapses. A synapse is a small gap across which one neuron will send information to another neuron in the form of electric current. The end of the neuron that is sending the signal is the axon, and the receiving end of the neuron is the dendrite. The act of sending a current across a synapse is referred to as neuron firing, and it occurs when the electric potential of the membrane, V m , reaches a certain threshold voltage, causing a spike in the voltage. This activity is caused by the moving of ions into and out of the neuron through the cell membrane, specifically sodium and potassium ions. At the threshold voltage, sodium gates open allowing sodium ions to flood in until a peak voltage is reached. It is at this time that the potassium gates open allowing potassium ions to flow out of the neuron, causing the potential to return to its resting state. Figure 1 shows a simple model of a neuron.
Controlling the firing rate of a neuron has numerous applications to various neurological diseases. It has already been shown that control systems can be used to help people with seizures [1][2][3][4][5][6] and to artificially recreate the experience of touch 7 . Another application is the control of neural oscillation and the blocking of neuron firing, which has been shown to be helpful in the treatment for diseases such as Parkinson disease or epilepsy [8][9][10] . Controlling a neuron to affect the motor functions of animals such as molluscs and insects is another application 11,12 . Additionally it has been shown that controlled electrical stimulation as a treatment for migraines 13 . In general, electrical stimulation to control action potential can potentially restore functions that an individual has lost, such as vision or hearing [14][15][16] .
Various researchers have considered the problem of controlling the dynamics of neurons in great detail. One method is to base the model around the Hodgkin-Huxley (HH) equations, which are four coupled and nonlinear differential equations that describe a change in the membrane voltage of a neuron [17][18][19][20][21][22] . Other researches have instead designed their models around the integrate-and-fire method [23][24][25] .
As to the controllers themselves, there has been work into using open loop and closed loop strategies to control how a neuron fires (for example, [26][27][28][29][30][31][32][33][34][35][36][37]. Some researchers have used model predictive control, a form of linear feedback control, to design their controller 17 , while others have used feedforward control 8,38 . Also used has been a series of adaptive controllers which allows one to tune the controller at different stages of design 38 .
For simplicity, this paper will be focusing on the control of only a single neuron's firing; the Hodgkin-Huxley equations model the behavior of a neuron and this firing. Then controllers will be designed using methods such as lead/lag and dynamic inversion.
The goal of this paper is to design and compare various model based controllers on the Hodgkin-Huxley model. These controllers and the response of the system are simulated and compared to determine which gives the most efficient response with the least error. This has not been done for the Hodgkin-Huxley model and thus is a novel approach.

Figure 1.
A simplified diagram of a neuron. As described, once the cell membrane reaches the threshold voltage, ions channels open allowing sodium and potassium ions to flow in and out, respectively. This results in the membrane voltage to spike, and an electric current moves across the synapse and is received by another neuron; this is how information is sent from neuron to neuron. The axon is the part of the neuron that carries the current away from the neuron, while the dendrite is what receives the current.

The Hodgkin-Huxley Model
As mentioned, the HH equations are four coupled non-linear differential equations which describe the change in V m . The equations take this V m to be dependent on the input current, I in j . Figure 2 shows the circuit representation of the neuron as described. The current of the flow of potassium ions is described by a conductance of G K , a voltage of E K which describes the potassium ion separation, and the term n which describes the behavior of the potassium gates opening. Similarly, we can define terms for sodium: G Na , E Na , m, and h. Notice that to describe the opening of sodium gates, two terms m and h are required. The terms G m and V leak describe the current of the remaining ions. The conductances depend on the voltage, while the current sources do not. A capacitor is included in parallel to represent the cell membrane, which behaves as a capacitor. Summing the currents in each component of the circuit gives Equation 1: with constants defined as G Na = 120 mS cm 2 , 613mV , and C m = 1 mF cm 2 . This equation is clearly non-linear, since m, n, and h are all variables representing the potassium gate activation, sodium gate activation, and sodium gate inactivation, respectively. They take on continuous values between zero and one representing the fraction of the channels in the corresponding state. They are described by first-order differential equations that depend on parameters α n , β n , α m , β m , α h , and β h who each in turn depend non-linearly on the membrane voltage, V . The equations are given as: 2/16 Figure 2. The circuit representation of the current flowing through a neuron membrane as described by the HH equations. The total current in the neuron is the sum of the current of each component of the circuit: the capacitor, the potassium, the sodium, a leakage current (G m ), and the injected current, I. The conductances G K , G Na , and G m represent the ion channels and the voltage sources E K , E Na , and E leak represent the equilibrium potentials.

Linearization
To begin analysis of this system, it will be linearized. We shall define the state space x ∈ R 4 to be The system is nonlinear, thus it may behave differently around different points and linearizing around one operating points may not be enough to describe the behavior of the entire system. The system will therefore be linearized at multiple operating points In general, the linearized state space of the system can be represented in the forṁ where the input, u, is the injected current I in j , and

3/16
and C = 1 0 0 0 (7) where . The response of each linearized system to a current step of I in j = 8nA is shown in Figure  3. Clearly the behavior does not change significantly enough to justify using multiple operating points. Furthermore, as will be shown later in this paper, designing a controller for the system linearized around one operating point meets the desired control objectives on the nonlinear system. Thus controllers will be designed for the system linearized about the point .
Clearly this linearized system is stable, as the poles are all in the left hand plane. The plot for V o = −70 mV in Figure 3 also shows this. Thus controllers can be designed for this system.

Feedback Control Strategies Lead Compensator
The objective is to control the membrane voltage by injecting current into the neuron, with the ultimate goal being to cause an activation spike by having the membrane voltage reach the threshold voltage of about −55 mV. From the transfer function in Equation 9 it can immediately be seen that the system is stable, as the poles are all in the left hand plane. A lead compensator is designed with general form G d (s) = k c s z +1 s p +1 , where z is the zero, p is the pole, and k c is the gain. To begin, the Bode plot of G(s) is examined, shown in Figure 4. The phase margin here is already high, so the design will keep it as such. So the desired phase margin is PM ≥ 107 degrees, with a bandwidth frequency of w bw = w c ≥ 10 radians/second. This will increase the speed of the transient response because the crossover frequency is being increased. Unlike a lag compensator, a lead compensator does not take the steady state error into account during the design of the controller.
To begin the design of the compensator, Figure 4 is examined and the phase of the transfer function at the desired crossover frequency, denoted as ∡G( jw c ), must be found. This is found to be −88.3 degrees. Then the required phase lead, φ m , can be found and the lead compensator pole and zero can be calculated as follows: Now these equations can be solved to get Then k c is chosen to ensure that the crossover frequency is at the desired value. This is done as follows: The Bode plot of the compensated system is shown in Figure 5. It can be seen that the compensated system meets specifications. Clearly the system meets the phase margin and crossover frequency specification.

Lag Compensator
A lag compensator is also designed, with the same specification of phase margin of ≥ 107 degrees and a steady state error of 1%. We shall call D(s) = s+z s+p . Figure 4 can be used to find the frequency at −180 + PM, where PM is the desired phase margin 107 degrees. Lag compensation will decrease the steady state error of the system, but unlike a lead compensator is does not improve the speed of the transient response. The frequency at this phase is about 0.985 radians/sec, and the magnitude at this phase is M = 0.9708. This magnitude can be used to calculate K, since K = 1 M = 1 0.9708 = 1.03. Then a lag compensator is designed as follows: where e is the desired steady state error. From there: We can choose a small z, such as z = 0.001, and then get p = 0.00003469. Finally the lag compensator is The Bode plot of the compensated system is shown in Figure 6. It can be seen that the compensated system meets the phase margin specification set. The error specification will be discussed later. Now the lag zero and pole are found using a similar process as the lag compensator.  Figure 7, it is clear that the lead/lag compensated system meets the requirements for phase margin. The error specification will be discussed later.

Direct Nonlinear Control
Another controller will now be designed based on dynamic inversion, a method of direct nonlinear control, which does not require linearization of the system. The nonlinear system can be defined aṡ and where u is the control input that must be found. From equations 1 and 2, we can define the values as follows and C = 1 0 0 0 The idea of dynamic inversion is to transform the coordinates such that the transformed system appears linear, and then a controller can be designed on this linear looking system.
Then u can be calculated as follows: The tracking error is defined as where y * is the threshold voltage, denoted henceforth as V th . Then select a fixed gain K such thaṫ Here, the gain K will be defined to be 20, as this provides acceptable response time and steady state error. Using the definitions in equations 24 to 29,Ė can be found and plugged into equation 30.
Equation 32 can now be solved for u.
Plugging in the values for f , g, and C gives a final value of u as A V th of −55 mV will be used here. This controller can then be applied directly to the system, shown in Figure 15 and discussed in the next section.

Results
All of the controllers designed were modeled in Matlab and Simulink in closed loop feedback as described by Figure 8. This allows for the membrane voltage to be compared to a reference voltage, which is the input for the closed loop system.
When simulating the lead, lag, and lead/lag controllers on the non-linear system, the system is allowed to reach a steady state of slightly below the threshold voltage instead, specifically −62 or −67 mV, depending on the system. Then the reference voltage is increased to pass the threshold voltage and when this new steady state is reached the controller is turned off, allowing the spike to occur. Random Gaussian white noise is also added to this input current for all simulations, as can be seen on the graphs, although on some plots it is harder to see than on others. 9/16 Figure 8. The block diagram of the controller on the nonlinear system (HH equations). Four controllers were used in this way: lead, lag, lead/lag, and dynamic inversion. Random noise is added when simulating the controller on the nonlinear system.
For the simulations of the lead, lag, and lead/lag controllers on the non-linear systems, the reference voltage will be a ramp input. This allows for better tracking of these systems, and for the membrane voltage to reach steady state before allowing it to spike. Also, using a step function for reference voltage would create a large error immediately, which would in turn result in a large spike in the current. Using a ramp for the reference voltage will allow for the current to increase smoothly, as can be seen in Figures 10, 12, and 14.
The tracking of the lead compensated system on the linearized system can be seen in Figure 9. Clearly steady state is reached in a reasonable amount of time and with a low steady state error. This error is calculated using the final value theorem: e ss = lim s→0 sError (35) which can be expressed as Equation 17, where KD(s) is the transfer function of the controller. A steady state error was not specified in lead compensator design, but nonetheless the error is 3.8%. The lead compensated system on the non-linear system can be seen in Figure 10; this system achieves a steady state error of 2.58%. As described before, once the reference voltage is increased and the membrane voltage reaches the threshold voltage, the controller is turned off allowing the voltage to spike.  Figure 10. Lead compensator G d used to track a reference voltage that ramps to −62 mV directly on the nonlinear system of Equations 1 and 2. Once steady state is reached, the reference voltage is changed to −55 mV, and once this threshold voltage is reached the controller is turned off, allowing a spike to occur. The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.
For the lag compensated system, the tracking of the reference voltage can be seen in Figure 11 on the linearized system. The lag compensated system eventually reaches a steady state with an error of 1.0% as calculated using Equation 17, and this meets the specification. Also, the lag compensated system takes longer to reach steady state than the lead compensated system. The lag compensator on the non-linear system can be seen in Figure 12, and once again the voltage spikes once the controller is turned off. This system has a steady state error of 1.49%. Figure 11. Lag compensator G g used to track a reference voltage of −55 mV on the linearized system G. The resulting closed loop transfer function is G g G 1+G g G . The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.

Figure 12.
Lag compensator G g used to track a reference voltage that ramps to −67 mV directly on the nonlinear system of Equations 1 and 2. Once steady state is reached, the reference voltage is changed to −55 mV, and once this threshold voltage is reached the controller is turned off, allowing a spike to occur. The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.
When the reference voltage is applied to the lead/lag compensated system, it can be seen through Equation 17 having the error be about 1.0%, which meets the specification. Also, the transient response is observed to be faster than that of lag. The tracking of this system can be seen in Figures 13 and 14 on the linearized and non-linear systems, respectively. The non-linear system achieves a steady state error of about 0.00% and spikes after the controller is turned off as described before.  Lead/lag compensator G lead/lag used to track a reference voltage that ramps to −62 mV directly on the nonlinear system of Equations 1 and 2. Once steady state is reached, the reference voltage is changed to −55 mV, and once this threshold voltage is reached the controller is turned off, allowing a spike to occur. The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.
Finally the tracking of the reference voltage for the dynamic inversion controller can be seen in Figure 15. The voltage reaches a steady state at the reference voltage with perfect tracking, and then like for all of the other systems, turning off the controller causes the voltage to spike. Figure 15. The dynamic inversion controller u used to track a reference voltage of −55 mV on the nonlinear system of Equations 1 and 2 (K=20). The controller is turned off after steady state is reached, allowing a spike to occur. The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.

Discussion
An activation spike in a neuron is when the voltage of the cell membrane reaches a certain value, causing current to move from the axon of that neuron to the dendrite of another across the synapse. In this way information is moved from neuron to neuron. This spike can be controlled, meaning that the spike can be triggered at will, a method that has a variety of purposes. As mentioned before, it can be used as a treatment for various neurological diseases, as well as helping to restore sight or hearing to those who are impaired. Four controllers were designed with the purpose of forcing a neuron to have this activation spike. Specifically the controllers are meant to determine the necessary current to be injected into the neuron cell membrane to reach the threshold voltage. Simulations of these controllers were done on a non-linear model of the neuron known as the Hodgkin-Huxley equations, as well as a linearized version of this system. These simulations were meant to show that the neuron could be forced to spike, and they did so.
As discussed, the controllers on the linear system all achieve tracking with low steady state error and a fast transient response, detailed in Table 1. The only exception is the lag compensator, which is the only controller with a slow transient response which is to be expected for this type of controller. On the nonlinear system, Table 2 shows that all of the controllers simulated were able to achieve low steady state error, and the dynamic inversion controller achieved a settling time comparable to that of the lead and lead/lag controllers on the linear system. However, for the first three controllers, the system was forced to track the ramp input for 20 milliseconds, and thus that is the time it takes those systems to reach steady state. The controllers on the nonlinear system have more error than the same controllers on the linearized system; this is because the controllers were designed on the linearized system, not the nonlinear system. However their error is still sufficiently low as to be acceptable.
With successful simulated control of a neuron's spiking, the next steps would be to study them on actual subjects, such as mice or insects. Figure 1 A simpli ed diagram of a neuron. As described, once the cell membrane reaches the threshold voltage, ions channels open allowing sodium and potassium ions to ow in and out, respectively. This results in the membrane voltage to spike, and an electric current moves across the synapse and is received by another neuron; this is how information is sent from neuron to neuron. The axon is the part of the neuron that carries the current away from the neuron, while the dendrite is what receives the current. The circuit representation of the current owing through a neuron membrane as described by the HH equations. The total current in the neuron is the sum of the current of each component of the circuit: the capacitor, the potassium, the sodium, a leakage current (Gm), and the injected current, I. The conductances GK, GNa, and Gm represent the ion channels and the voltage sources EK, ENa, and Eleak represent the equilibrium potentials. The system was linearized about four different operating points to determine how the behavior of the system changed. The four points were chosen between the resting voltage -70 mV and the threshold voltage -55 mV. From this four transfer functions were calculated, and their response to a step current (Iin j = 8nA) is plotted here.

Figure 4
The Bode plot of the transfer function of the linearized system, G, is shown here. The Bode plot was used in designing lead, lag, and lead/lag controllers.

Figure 5
The Bode plot of the lead compensator on the linearized system. The resulting open loop transfer function is GdG. Clearly the system meets the phase margin and crossover frequency speci cation.

Figure 6
The Bode plot of the lag compensator on the linearized system. The resulting open loop transfer function is GgG. Clearly the system meets the phase margin speci cation.

Figure 7
The Bode plot of the lead/lag compensator on the linearized system. The resulting open loop transfer function is GleadlagG. Clearly the system meets the phase margin and crossover frequency speci cation.

Figure 8
The block diagram of the controller on the nonlinear system (HH equations). Four controllers were used in this way: lead, lag, lead/lag, and dynamic inversion. Random noise is added when simulating the controller on the nonlinear system. Please see the Manuscript PDF le for the complete gure caption.

Figure 10
Lead compensator Gd used to track a reference voltage that ramps to -62 mV directly on the nonlinear system of Equations 1 and 2. Once steady state is reached, the reference voltage is changed to -55 mV, and once this threshold voltage is reached the controller is turned off, allowing a spike to occur. The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.

Figure 11
Please see the Manuscript PDF le for the complete gure caption.

Figure 12
Lag compensator Gg used to track a reference voltage that ramps to -67 mV directly on the nonlinear system of Equations 1 and 2. Once steady state is reached, the reference voltage is changed to -55 mV, and once this threshold voltage is reached the controller is turned off, allowing a spike to occur. The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.

Figure 13
Please see the Manuscript PDF le for the complete gure caption.

Figure 14
Lead/lag compensator Glead/lag used to track a reference voltage that ramps to -62 mV directly on the nonlinear system of Equations 1 and 2. Once steady state is reached, the reference voltage is changed to -55 mV, and once this threshold voltage is reached the controller is turned off, allowing a spike to occur.
The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.

Figure 15
The dynamic inversion controller u used to track a reference voltage of -55 mV on the nonlinear system of Equations 1 and 2 (K=20). The controller is turned off after steady state is reached, allowing a spike to occur. The injected current produced by the controller with the addition of random Gaussian white noise is also plotted.