Search tips
Search criteria 


Logo of springeropenLink to Publisher's site
Journal of Computational Neuroscience
J Comput Neurosci. 2010 June; 28(3): 375–387.
Published online 2010 February 5. doi:  10.1007/s10827-010-0215-x
PMCID: PMC2880706

A model of feedback control for the charge-balanced suppression of epileptic seizures


Here we present several refinements to a model of feedback control for the suppression of epileptic seizures. We utilize a stochastic partial differential equation (SPDE) model of the human cortex. First, we verify the strong convergence of numerical solutions to this model, paying special attention to the sharp spatial changes that occur at electrode edges. This allows us to choose appropriate step sizes for our simulations; because the spatial step size must be small relative to the size of an electrode in order to resolve its electrical behavior, we are able to include a more detailed electrode profile in the simulation. Then, based on evidence that the mean soma potential is not the variable most closely related to the measurement of a cortical surface electrode, we develop a new model for this. The model is based on the currents flowing in the cortex and is used for a simulation of feedback control. The simulation utilizes a new control algorithm incorporating the total integral of the applied electrical potential. Not only does this succeed in suppressing the seizure-like oscillations, but it guarantees that the applied signal will be charge-balanced and therefore unlikely to cause cortical damage.

Keywords: Epilepsy, Seizure, Feedback, Control, Electrocorticogram, Human


The recurrent, unprovoked seizures associated with epilepsy can have a devastating effect on those with this disorder. Basic parts of every day life such as driving and obtaining employment become very difficult. While many people with epilepsy can control their seizures with medication, roughly thirty percent do not respond to this type of treatment and therefore seek out alternatives such as surgery (The Epilepsy Foundation 2009). The surgical procedure involves resecting the seizing portion of the brain, often part of the cortex, while avoiding any areas that provide vital functions such as speech, memory, and vision (The Epilepsy Foundation 2009). Because this is a very invasive procedure that does not guarantee success, other alternatives are being investigated. One of these is automatic feedback control, where subdural electrodes on the cortical surface would detect the seizure and apply an electrical signal to disrupt the abnormal electrocorticogram (ECoG) activity.

This method of treatment is currently being studied by experimentalists. It has been shown that the application of electric fields to rat cortex in vitro can modulate the behavior of seizure-like waves (Richardson et al. 2005). In vivo experiments on rats demonstrated that stimulation via proportional feedback can temporarily suppress seizure activity (Gluckman et al. 2001). A subsequent set of experiments showed that an increase in the amplitude of the proportional control feedback gain corresponds to decreases in both seizure amplitude (measured as a reduction in the amplitude variance) and Teager energy (Colpan et al. 2007). The Teager energy can be reduced by a decrease in the amplitude of a signal or a lowering of the frequencies in its power spectrum.

While these experiments suggest further exploration, they are all restricted to animal implementation. Therefore, to gain insight into the feasibility of human use, we turn first to mathematical models. The stochastic partial differential equation (SPDE) model of the cortex used here can support seizure-like oscillations that are qualitatively and quantitatively similar in frequency of maximum power and propagation speed to those seen in humans with epilepsy (Kramer et al. 2005). We previously demonstrated that various methods of feedback control can suppress these simulated seizures (Kramer et al. 2006), and we added the capability of looking at spatial properties of feedback control, such as electrode size and spacing (Lopour and Szeri 2008).

Our aim in the present work is to make this model more biologically relevant by refining the representation of feedback control. This will facilitate future comparison with experimental data. There are four key improvements in our approach that are presented in this article:

  1. We verify the strong convergence of the numerical solution to the SPDE model.
  2. Based on the studies of convergence, we utilize a smaller step size in our simulations, thereby allowing the inclusion of a more detailed electrode profile.
  3. We develop a better motivated model of the signal measured by an electrode on the cortical surface. This model is used to calculate the applied electrical signal for feedback control.
  4. Feedback control is performed with a new algorithm incorporating an integral component. This ensures that the applied signal is charge-balanced, which is thought to minimize damage of cortical tissue.

The organization of the paper is as follows. We first briefly review the SPDE cortical model that will be used in our simulations (Section 2), and then we discuss the convergence of its numerical solutions (Section 3). Next, we present the new model for electrode measurements (Section 4). Finally, we incorporate these improvements into simulations of feedback control, while implementing a new integral control law (Section 5).

Cortical model

To model the electrical activity of the human cortex, including seizure waves, we choose a set of stochastic PDEs that has been developed and adapted over the past decade (Liley et al. 1999, 2002; Steyn-Ross et al. 2003). The mesoscale nature of this model makes it well-suited to EEG-based applications such as epilepsy, sleep (Wilson et al. 2006), and anesthesia (Steyn-Ross et al. 2004) because it is based on length scales similar to commercial electrode arrays. It is a mean-field model, meaning that all of its variables represent spatially averaged properties of populations of neurons. This is similar to the manner in which an electrode provides a measurement based on the collective behavior of many neurons.

In 2006, the equations were restated in a dimensionless form by Kramer et al. (2007). This is the formulation of the model we will use here; it is a system of eight coupled nonlinear PDEs with stochastic inputs:

equation M1
equation M2
equation M3
equation M4
equation M5
equation M6
equation M7
equation M8

All variables are dimensionless and are functions of time (equation M9) and one spatial dimension (equation M10). The equation M11 state variable is the mean soma potential for a neuronal population, while equation M12 represents postsynaptic activation due to local, long-range, and subcortical inputs, and equation M13 is a long-range (corticocortical) input. The subscripts e and i denote affiliation with the excitatory and inhibitory neuron populations, respectively; variables with two subscripts represent the transmission of information from one population to another, e.g. equation M14 is the postsynaptic activation of the excitatory population due to inputs from the inhibitory population. In Eq. (1) we have added the variable u to represent the signal applied by a cortical surface electrode for feedback control. This will be discussed further in Section 5. For descriptions of all model variables and parameters, please refer to Table 1.

Table 1
Dimensionless variables and parameters of the SPDE cortical model

To appreciate the model as a whole, let us first look at the equations governing the excitatory neuron population, depicted graphically in Fig. 1. Equation (1) for equation M42 is reminiscent of the leaky integrate-and-fire model of a neuron, where the derivative of the membrane potential equals the resting potential minus the membrane potential plus any existing current inputs (Dayan and Abbott 2001). Here, the resting potential is “1” due to the dimensionless nature of the system. The inputs are equation M43 and equation M44, which evolve according to (3) and (5), respectively, based on three types of synaptic input: local, long-range, and subcortical.

  1. Local inputs, such as those from within the same macrocolumn, are represented by terms of the form equation M45, where equation M46 is a dimensionless sigmoid function:
    equation M47
    This converts the mean soma potential of the excitatory population to its mean firing rate.
  2. Long-range inputs represent signals from other cortical macrocolumns and are defined by equation M48. The behavior of this variable is governed by (7). Note the similarity of this equation to the standard PDE wave equation; the idea that cortical tissue can support wave propagation is central to our simulation of epileptic seizures.
  3. Subcortical inputs are predominantly from the thalamus and contain both constant (Pee) and stochastic (equation M49) parameters. We define the stochastic term by
    equation M50
    where αee is a constant and ξ1 is zero mean, Gaussian white noise in time and one spatial dimension. When the SPDEs are solved numerically, the cumulative effect of this stochastic process will be Brownian motion. To ensure that the properties of this signal remain constant regardless of step size, we scale the discrete randomly generated numbers R(m,n) by the simulation time step:
    equation M51
    The variables m and n are indices of space and time, so a single point is represented by spatial position equation M52 at time equation M53. Note that the discrete update form of (11) is
    equation M54
    which will be used in all numerical experiments.
Fig. 1
Flow chart representation of model Eqs. (1), (3), (5), and (7), which govern the excitatory population. The boxes describe the physiological significance of the model variables and parameters listed beneath them. Note that feedback occurs through ...

Thus, Eqs. (1), (3), (5), and (7) govern the excitatory population, while the remaining equations represent the inhibitory population and have exactly the same form. Together they compose the full cortical model.

There are several parameters that are especially relevant to the following numerical studies. The parameter Pee represents input from the population of subcortical excitatory neurons (such as those in the thalamus), and Γe denotes the influence of synaptic input on the mean soma potential. Changes in these parameters allow for transition between normal cortical function (Pee = 11.0 and equation M58) and the hyperexcited “seizure” state of the SPDE model (say, Pee = 548.0 and equation M59). At low levels of excitation corresponding to low levels of Pee, the mean soma potential of the excitatory neurons equation M60 produces random fluctuations similar to those seen in an EEG measurement. However, at increased levels of subcortical excitation, the simulated cortex develops large amplitude seizure-like oscillations. Our goal is to suppress this pathological behavior via feedback control consisting of measurements from the cortical surface and the application of a potential based on those measurements.

Strong convergence of numerical solutions

Before performing simulations of feedback control, we must ensure that we can obtain accurate numerical solutions to this system of SPDEs. We will use a predictor-corrector algorithm written in MATLAB, so the accuracy of the solution will be determined by our choice of step sizes in space and time. In addition to considering the system equations and solution method, we shall assume that a typical cortical surface electrode is of order 1 cm in diameter. While the previously used step size of 7 mm (Lopour and Szeri 2008; Kramer et al. 2006) may have accurately solved the differential equations, it was not small enough for sufficient spatial resolution of the behavior of the cortical tissue underneath the electrodes. We will use much smaller step sizes in order to achieve both of these objectives.


To determine the magnitude of these step sizes, we will examine the strong convergence of solutions to the cortical model. This will be accomplished by generating equivalent Brownian paths at several step sizes and demonstrating that the solutions converge as the step size decreases (note that this differs from weak convergence, which looks at the expected value of the solution over all possible Brownian paths) (Higham 2001). This task will be complicated by the fact that both the stochastic inputs and the numerical solutions vary in space and time.

Recall from (12) that the grid of stochastic inputs is defined by R, which consists of M independent Brownian paths, each of length N. This corresponds to M points in space at a reference step size of Δx = Δx0 and N points in time at a reference step size of Δt = Δt0. Therefore, we denote individual points by R(m,n), ∀ m = 1,...,M and ∀ n = 1,...,N. Then we can represent the same Brownian paths at a coarser step size 2Δt0 by adding together every two adjacent elements in time:

equation M61

where n = 1,...,N/2 (Gaines 1995). We do not need any special scaling factors here because this combination of neighboring terms is consistent with the definition of a Brownian path. Similarly, we can represent the stochastic input at step size 2Δx0 by adding together adjacent elements in space and scaling to keep the variance constant (Gaines and Lyons 1997):

equation M62

where m = 1,...,M/2. The factor equation M63 is necessary because the stochastic inputs are independent in the spatial direction; it will be used for the relative scaling of the inputs at different step sizes for the purpose of determining convergence, but will not be present in a typical simulation of feedback control.

Now we can directly compare numerical solutions at decreasing step sizes (e.g. Δx = x0, 2Δx0, Δx0) under equivalent stochastic inputs. We want to verify that the solution converges as we approach Δx0.


First, we look at the convergence in time using the method described above. We remove the spatial terms from the cortical model to reduce it to an ODE and then perform simulations with decreasing values of the time step. These indicate that the solution converges around Δt = 5 ×10  4 s (Fig. 2). The two smallest time steps in the figure, Δt = 5 ×10  4 and 2.5×10  4 s give very similar results for he.

Fig. 2
Convergence of numerical solutions as the time step is decreased using the method described in Section 3.1. Here the spatial terms in the model have been removed to reduce it to an ODE. We used equation M64 s, and we plotted the solutions for Δ ...

Next, we study the convergence in space. We begin with a Δx that is smaller than 1 cm because we desire to resolve the solution across an electrode. As Δx decreases, we see that the accuracy of the solution improves; however, this does not give us a clear indication of which step size to choose. The amount of improvement seems to be the same for each reduction in Δx. We solve this problem by looking at the numerical worst-case scenario—a sharp transition between uncontrolled cortex and a single electrode with proportional feedback. We then choose Δx based on its ability to resolve this sharp spatial change (Fig. 3). While this figure shows that the differences between the step sizes are still subtle, it seems that the largest one, Δx = 0.448 mm, does not provide enough detail to show the sharp transition between cortex and electrode. The smaller step sizes appear to be more accurate and provide very similar solutions. Because it will be less computationally intensive to use Δx = 0.224 mm (or equation M66 dimensionless), we choose this as the step size for our simulations.

Fig. 3
Convergence of numerical solutions over an electrode nonlinearity as the spatial step size is decreased using the method from Section 3.1. The full simulation spanned 22.4 mm, but here we show only 3 mm of uncontrolled cortex and ...

In order to verify this in the typical case with no feedback control, we plot numerical solutions with Δx = 0.448, 0.224, and 0.112 mm (Fig. 4). Because the spatial solution appears to have converged at these step sizes, our choice of Δx = 0.224 mm is valid. The last task is to choose a value of Δt. We would like to use the largest time step for which the solution converges because this will result in the shortest computation time; this value is equation M70 (dimensionless) or Δt = 4×10  6 s. Because we already showed that accurate solutions can be obtained with much larger values of Δt, this is an acceptable choice.

Fig. 4
Convergence of numerical solutions as the spatial step size is decreased using the method described in Section 3.1. Here we used Δx0 = 0.112 mm, and we plotted the solutions for Δx = ...

Therefore, the step sizes used in all of the following cortical simulations will be Δx = 0.224 mm and Δt = 4×10  6 s.

Model of feedback control

Our previous simulations of feedback control (Kramer et al. 2006; Lopour and Szeri 2008) utilized two key assumptions: 1) the signal measured by an electrode on the cortical surface is proportional to he, the mean soma potential of the excitatory neuron population, and 2) a voltage applied to the surface of the cortex via electrode directly affects the average soma voltage in that region. The first assumption allows us to define the control effort u in terms of equation M73 (in this case, equation M74 would represent the measured voltage), and the second assumption implies that the expression for u can be added directly to the SPDE model in Eq. (1). While the latter assumption appears to be valid, there is evidence that we cannot write u as an explicit function of equation M75 as the first assumption suggests. It is likely that the voltage sensed by a surface electrode is different than the averaged soma voltage, equation M76.

First, it is important to realize that the signal measured by an electrode is a function of the extracellular currents in the tissue, rather than the intracellular somatic potential (Nunez and Srinivasan 2006). We define the signal sensed at a point on the cortical surface to be equation M77. Then, to understand the difference between equation M78 and equation M79, we consider a pyramidal neuron in the cortex with one excitatory synapse as shown in Fig. 5. Say that the pyramidal neuron receives excitatory input due to a proximal synapse in layer 4 (Fig. 5(a)); this will cause intracellular flow of ions that will induce a current dipole with sources (+) on the apical dendrite near the surface and sinks (−) near the soma. The surface electrode equation M80 will sense the extracellular current source near the surface and will thus depolarize. The soma potential equation M81 will also depolarize due to the excitatory input; therefore, in this case, both equation M82 and the surface electrode show a depolarization. On the other hand, suppose that the pyramidal neuron receives excitatory input due to a distal synapse in layer 1 (Fig. 5(b)). Because the input is still excitatory, the neuron will depolarize, and this will be reflected in the soma potential equation M83. However, this input will cause an extracellular current dipole with reverse polarity; there will be a source (+) near the soma and current sinks (−) near the surface. This means that the voltage sensed by the surface electrode will show a hyperpolarization. Therefore, in this case, the deflection of equation M84 and the signal seen by the surface electrode are different (Kandel et al. 2000).

Fig. 5
Relationship of the sensed signal equation M85 to the mean soma potential equation M86. In (a), an excitatory input synapses near the soma, causing a depolarization in equation M87. Similarly, the dipole current generated by this synaptic event involves current sources near the cortical ...

This implies that we should no longer use equation M91 as the measured electrical potential in our expression for the control effort u. Instead, the measurement will be a function of the currents in the cortex due to synaptic inputs, denoted as equation M92. We will refer to this as the sensed signal. Then, for the purposes of feedback control, the applied effort u will be a function of equation M93.

Basic form of model

To determine the composition of the sensed signal, we must consider the extracellular current flows due to three types of synaptic input (local intracortical input, long-range corticocortical input, and subcortical input). To do this, we need to know whether the inputs are excitatory or inhibitory and whether they synapse near the surface or near the soma. In what follows, we take care to distinguish depolarization in the electrode measurement from depolarization of the soma.

  1. Local intracortical inputs, equation M94 and equation M95. Within a cortical macrocolumn, excitatory synapses tend to occur close to the surface, while inhibitory synapses are located near the soma (Nunez and Srinivasan 2006; Spruston 2008). Thus, the excitatory inputs are as depicted in Fig. 5(b), and the inhibitory inputs have the geometry of Fig. 5(a) but with the opposite sign (because Fig. 5 depicts excitatory inputs). Each of these configurations will cause a hyperpolarization in the electrode measurement; therefore, both terms will have negative signs in the measurement model: equation M96, where A and B are positive constant weights to be determined.
  2. Long-range corticocortical input, equation M97. Corticocortical inputs are exclusively excitatory (Steyn-Ross et al. 2003; Liley et al. 2002; Nunez and Srinivasan 2006) and tend to synapse near the surface (Nunez and Srinivasan 2006; Spruston 2008). More specifically, layers 2 and 3 of the cortex seem to have a higher density of corticocortical inputs (Kandel et al. 2000; Nieuwenhuys 1994). As shown in Fig. 5(b), this input type will cause a hyperpolarization in the electrode signal and will thus have a negative sign in the measurement model: equation M98, where C is a positive constant weighting factor to be determined.
  3. Subcortical inputs, equation M99 and equation M100. While the distribution of these synapses is not clear-cut, it seems that subcortical inputs terminate most densely in layer 4 near the soma (Kandel et al. 2000; Nieuwenhuys 1994). Because Pee is an excitatory input of the type shown in Fig. 5(a), it will have a depolarizing effect on the electrode measurement; therefore, we give it a positive sign: equation M101, where D is a positive constant weight. On the other hand, Pie is inhibitory and will thus have a hyperpolarizing effect: equation M102, where E is a constant weighting factor. The values of D and E are to be determined.

Incorporating all three types of synaptic input gives us this basic expression:

equation M103

with A, B, C, D, E of positive sign but (so far) unknown magnitude. The consequence of these inputs is only evident after synaptic transmission. Therefore, we include a rate constant for this process by using an equation similar to that of equation M104 in the SPDE model. Let equation M105 represent the current measured at the cortical surface and Tm represent a rate constant. Then

equation M106

where A, B, C, D, E, and F are positive constant weights. We choose Tm = 12.0 to match the rate constant of the excitatory population Te. The values of A through E will depend on the number of synapses of each type and the average distance of the synapse from the soma. The coefficient F is a gain parameter that will scale the magnitude of all the synaptic inputs; this ensures that they have the appropriate amount of influence over the electrode measurement equation M107. In addition, we can think of F as containing the effective resistance of the cortex. Recall that the electrode measurement is determined by currents in the cortex, yet the components on the right side of (16) are based on voltages. Because the currents produced by these voltages can be calculated with Ohm’s Law (Kandel et al. 2000), the gain parameter F provides the necessary conversion.

To complete the model of the electrode measurement equation M108, we must account for the reversal potential of the cortical neurons. This determines the direction of current flow associated with the inputs described above (we previously assumed that the neurons were at resting potential). We once again take our cue from the SPDE model and define

equation M109

Thus, (16) and (17) comprise a complete model of the potential sensed by a cortical surface electrode, equation M110. In our simulations of feedback control, the applied electric field u will be a function of this variable.

Estimation of coefficients

We have not yet assigned numerical values to the coefficients A through E. To do this, we think of them as the percentage of pyramidal neuron synapses due to each source. For example, A will represent the percentage of synapses on any given pyramidal neuron that come from other excitatory neurons in the same macrocolumn. There are three physiological relationships that allow us to determine these values:

  1. The number of synapses on pyramidal cells due to local cortical neurons is roughly equal to the number of synapses due to cortical neurons in other macrocolumns or in the contralateral hemisphere (Abeles 1991). This implies that A + B = C.
  2. Approximately 98 percent of synapses on pyramidal cells are corticocortical, while 2 percent are thalamocortical (Abeles 1991; Nunez and Srinivasan 2006). This implies that A + B + C = .98 and D + E = .02.
  3. Roughly 90 percent of all cortical synapses are excitatory and 10 percent are inhibitory (Abeles 1991; Braitenberg and Schüz 1998). This implies that A = 9B and D = 9E.

After solving the above equations, we account for the fact that synapses near the soma will have a greater influence on the electrode measurement (Nunez and Srinivasan 2006) by multiplying B, D, and E by a factor of two.

While this method of estimation may seem crude, others have achieved similar results through more detailed probabilistic analysis (Liley and Wright 1994). Both are listed in Table 2 for comparison; note that each set of coefficients has been scaled to add to 1. In simulation, when comparing the two sets of coefficients A through E, the only difference in the resulting equation M111 appears to be an offset. Because offsets are not reflected in EEG measurements, this difference is inconsequential. Hence, two completely independent approaches provide essentially equivalent coefficients for the sensed signal.

Table 2
Values of the coefficients for Eq. (16)

Verification of full model

We can verify our model by returning to the hypothetical pyramidal neuron described at the beginning of Section 4 and based on Box 46-1 in (Kandel et al. 2000). As before, say that we are modeling the electrode measurement of a pyramidal neuron with excitatory inputs in both cortical layers 1 and 4. If a majority of the inputs occur in layer 1 near the cortical surface as in Fig. 5(b), equation M112 and equation M113 will have similar dynamics, but a hyperpolarization in one signal will be a depolarization in the other; the signals will be negatively correlated. This behavior is seen in our model for equation M114. If we run a simulation with the typical set of parameters (where local and corticocortical connections dominate because Pee is low), we see that equation M115 and equation M116 are negatively correlated (top of Fig. 6). On the other hand, when the strongest input is near the soma in layer 4 as in Fig. 5(a), it will have the same effect on both signals, and they will be positively correlated. This, too, is demonstrated by the measurement model. We can simulate a large excitatory input near the soma by increasing the value of Pee; when we run the simulation with this change, we see that equation M117 and equation M118 become positively correlated (bottom of Fig. 6). Thus, the measurement model accurately reproduces the physiological effects of varying cortical inputs, and we will use it in subsequent simulations of feedback control. The applied electric field u will be a function of equation M119 as opposed to equation M120.

Fig. 6
Comparison of equation M121 (solid) and equation M122 (dashed) with two sets of parameters. At typical levels of excitation where Pee = 11.0 and equation M123, the two signals are negatively correlated (top subfigure), as predicted by our physiological analysis. This negative ...

For reference, we compare equation M127 and equation M128 at seizure parameters (equation M129 and Pee = 548.0) in Fig. 7. In this case, the signals have a large negative correlation, but very similar dynamics. This suggests that it will be possible to suppress seizures with feedback control based on equation M130 (as it was with equation M131 in (Kramer et al. 2006)), although we may need to use a gain of the opposite sign.

Fig. 7
Comparison of the electrode measurement hm (mV, solid) and the mean soma potential he (mV, dashed) at levels of subcortical excitation that cause seizure-like oscillations, Pee = 548.0 and equation M132. Here we see that the two signals are negatively ...

Simulation of integral control

In choosing a function to represent the applied electric field u, we start with the concept of proportional feedback control. This is the simplest and most common type of control—intuitively, the applied effort should be proportional to the error between the measured signal and its desired value. In this case, we can define proportional control as

equation M135

where a(x,t) is the control gain and equation M136 is the measured electrode potential; it is calculated by taking the spatial average of equation M137 (17) under each electrode. The parameter b is a constant offset that can be tuned to achieve the desired equilibrium value of equation M138. While this control algorithm is able to suppress the seizure-like oscillations of the model, it would be difficult to implement safely. Whenever stimulation is applied, it is important that the process be chemically reversible in order to prevent damage due to the production of new chemical species. There is a threshold for reversibility called the “reversible charge injection limit,” which signifies the maximum allowable charge injection before the polarity is reversed (Robblee and Rose 1990). Because a proportional controller does not penalize the amount of effort used (i.e. the magnitude of the applied electric field), it relies on large signals of only one sign, which would exceed this threshold over time. The chemical processes associated with this type of stimulation would therefore be irreversible and damaging to cortical tissue. The simulation results based on this type of control have been presented in previous publications (Kramer et al. 2006), and we do not repeat them here.

To improve on this method, we may consider adding a derivative or integral component to the controller, or even using all three terms to create a proportional-integral-differential (PID) controller (Franklin et al. 2002). The derivative term increases or decreases the control effort based on the rate of change of the error. This can reduce the response time of the controller because the derivative term “anticipates” the behavior of the system. Simulation results with a PD controller were presented in Kramer et al. (2006).

Because the differential controller utilizes the same harmful voltages mentioned previously, we choose to implement a controller with an integral term:

equation M139

Here, c(x,t) is another gain term. It will be negative, meaning that this new term will oppose the total integral of the applied voltage u. If the integral of u is positive, it will add a negative component to the applied voltage, and vice versa; in this way, it pushes the integral of u to zero. In other words, it forces the applied signal to be charge-balanced and thus safe for cortical tissue. Because we have also included the proportional control term, this feedback setup will still suppress the seizure oscillations. Note that this is different than traditional integral control, which is based on the integral of the error between the desired value of the signal and its actual value.

In addition to adding the integral term and using the new measurement model equation M140, we incorporate the step sizes determined in Section 3. Because we are using smaller increments of space, we can make one further improvement: we add a more detailed electrode profile to the feedback simulation. We previously assumed that the electrodes maintained a constant profile across their surface while measuring or applying the stimulus (i.e. every point on the electrode sensed or provided the same value) and that no tissue beyond the edge of the electrode was affected by its activity. However it has been shown experimentally that this is not the case (Suesserman et al. 1991). Here, we take the first step towards a realistic electrode model by including a smooth falloff at the electrode edges. The falloff is incorporated into the function a(x,t), so that the gain varies between zero (over uncontrolled cortex) and a max (under the electrodes) via a hyperbolic tangent function. This gain function is used in the application of control to indicate that the influence of the electrode decreases with distance, and it is also applied during sensing with a max  = 1 to indicate that the cortical tissue has less impact on the electrode measurement as distance increases. We have not yet included any variation over the surface of the electrode, but this is certainly an adjustment that can be considered in the future.

With this approach, we are in a position to simulate the suppression of seizure waves using total integral feedback control as defined in (19). The results shown here were generated using an Intel Core 2 Duo 2.13 GHz processer, and the calculations took roughly 5 min at the highest spatial and temporal resolutions. The code was written and executed in MATLAB and has been provided as supplementary online material (Online Resource 1).

Figure 8 shows an example of the seizure-like oscillations, represented by the mean soma potential of the excitatory population he as it varies in space and time. More specifically, we have simulated an uncontrolled strip of cortex 200 mm long over 0.5 s, and the value of he [mV] is represented by grayscale. The seizure waves occur due to our choice of equation M141 and a Gaussian Pee distribution with a maximum value of 548.0. They spontaneously arise in locations with sufficiently high Pee (here, this “hot spot” is at x = 100 mm) and travel outward until the level of excitation is too low to support them. This is why the waves terminate before they reach the edges of the simulation space.

Fig. 8
(a) Seizure waves traveling on the simulated cortex, with parameters Pee = 548.0 and equation M142. The characteristics of the wave are determined by the distribution of Pee; here it is a Gaussian curve, so the wave starts in the center where Pee is ...

Figure 9 shows the effect of total integral feedback control on this seizure-like behavior. In the first half of the time interval, we see a seizure wave emanating from the center of the space; then, at 0.25 s the control is switched on, and the pathological seizure activity quickly disappears. In this case, the feedback stimulation is applied via five electrodes that are 11.2 mm across. The geometry of these electrodes can be clearly seen in Fig. 10, which shows the value of the applied signal u [mV] for the same simulation. Here, we see that the applied voltage is zero for the first half of the time interval, and then each electrode varies between positive and negative voltages once the controller is turned on. While the largest potential applied by a single electrode is roughly 60 mV, the total signal applied by each electrode is very close to zero, due to our choice of control law u. This is the desired result because it indicates that each electrode applies a balanced signal. Lastly, note that it is not necessary to have electrodes covering the entire length of the seizure wave. When we place the electrodes at the center of the “hot spot,” they are able to halt the outward motion of the wave. For more information on the parameters used to create the figures, please refer to the figure captions.

Fig. 9
(a) Plot of he [mV] in space and time. In the first 0.25 s, we see a seizure wave traveling on the simulated cortex because we have set the excitation parameters to Pee = 548.0 and equation M143. At 0.25 s (dashed line), the integral ...
Fig. 10
(a) Plot of applied effort u [mV] in space and time, corresponding to the simulation of feedback control in Fig. 9. This shows that the potential applied by the electrodes is zero until t = 0.25 s, when the controller is ...

Thus, in Figs. 9 and and1010 we have demonstrated that the new model for electrode measurements equation M145 can be used to suppress seizure waves via feedback control. If the controller u contains the integral of the total applied effort, then this can be done in a manner that is thought to be safe for cortical tissue.


Here we have presented several novel approaches to exploration of a model of feedback control for epileptic seizures in humans. We first verified the strong convergence of numerical solutions to the model of the cortex, paying special attention to discontinuities that may occur at electrode edges. This allowed us to choose appropriate step sizes for our simulations; because the spatial step size Δx was small relative to the size of the electrode, we were able to incorporate a more detailed electrode profile into the simulation. Then, based on evidence that the mean soma potential equation M146 cannot be used as the measurement for feedback control, we developed a new model equation M147 to represent the measurement of cortical surface electrodes. This model was based on the currents flowing in the cortex and was used for all simulations of feedback control. Those simulations utilized a new control algorithm containing the total integral of the applied potential u. Not only did this succeed in suppressing the seizure-like oscillations, but it guaranteed that the applied signal would be charge-balanced and therefore safe for cortical tissue.

Of course, there are always improvements to be made. In this work, we have assumed that each electrode can be simultaneously sensing and applying the control signal. This is not realistic; ideally, we would model separate electrodes for these two tasks, with the geometric properties chosen to match existing experimental setups. Also, as mentioned previously, it would be possible for us to improve the electrode profile used in simulation. Here, we incorporated a simple profile to demonstrate our capability to do so, but it would be more accurate to base our choice on existing theories of the potential difference across an electrode surface (Rubinstein et al. 1987). It may also be possible to account for electrochemical changes that occur in the vicinity. Lastly, we note that our use of the phrase “charge-balanced” should be taken lightly. Our controller measures and applies a voltage, so the integral term pushes the total voltage (over time) towards zero. Although, in concept, this is similar to having a charge-balanced signal, it does not guarantee that the applied charges will be balanced and safe. This could be remedied by utilizing a controller that measures cortical potential and applies a current. Not only would this be more accurate, but it would facilitate future comparisons with experiments, most of which are done in this manner (Colpan et al. 2007; Sunderam et al. 2006).

Validation via experimentation is only one of the many possible future directions of this work. For example, we could extend our model to two dimensions and use it to study seizure waves, or we could simulate experimental phenomena such as irregular, spiral, and plane cortical waves (Schiff et al. 2007). Related theoretical work has suggested that pre-processing of data using a Kalman filter can provide greater flexibility in the control of waves while minimizing the amount of energy needed to do so (Schiff and Sauer 2008). This concept could be readily applied to the simulations discussed here. Another possible avenue of future work is the investigation of spatial properties of our feedback model. We have the capability to do simulations with any number of electrodes at any size and spacing, which is a luxury not afforded to experimentalists. Theoretical work in this area may provide insight into important questions such as: how does the size of a seizure relate to the number and sizes of the electrodes needed to control it successfully? Where should the electrodes be placed for maximum effectiveness? What are the necessary resolutions for sensing and actuation? It is our hope that this work will act as a stepping stone to such intriguing questions.


This material is based upon work supported under a National Science Foundation Graduate Research Fellowship. We also extend special thanks to Bruce Gluckman and Steven Schiff for their time and valuable comments.

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Contributor Information

Beth A. Lopour, ude.yelekreb@ruopolhteb.

Andrew J. Szeri, ude.yelekreb.em@irezsa.


  • Abeles M. Corticonics: Neural circuits of the cerebral cortex. Cambridge: Cambridge University Press; 1991.
  • Braitenberg V, Schüz A. Cortex: Statistics and geometry of neuronal connectivity. 2. Berlin: Springer; 1998.
  • Colpan ME, Li Y, Dwyer J, Mogul DJ. Proportional feedback stimulation for seizure control in rats. Epilepsia. 2007;48(8):1594–1603. doi: 10.1111/j.1528-1167.2007.01073.x. [PubMed] [Cross Ref]
  • Dayan P, Abbott LF. Theoretical neuroscience: Computational and mathematical modeling of neural systems. Cambridge: MIT; 2001.
  • Franklin GF, Powell JD, Emami-Naeini A. Feedback control of dynamics systems. 4. Upper Saddle River: Prentice Hall; 2002.
  • Gaines JG. Numerical experiments with S(P)DE’s, chapter in Stochastic Partial Differential Equations. Cambridge: Cambridge University Press; 1995. pp. 55–71.
  • Gaines JG, Lyons TJ. Variable step size control in the numerical solution of stochastic differential equations. SIAM Journal on Applied Mathematics. 1997;57(5):1455–1484. doi: 10.1137/S0036139995286515. [Cross Ref]
  • Gluckman BJ, Nguyen H, Weinstein SL, Schiff SJ. Adaptive electric field control of epileptic seizures. Journal of Neuroscience. 2001;21(2):590–600. [PubMed]
  • Higham DJ. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review. 2001;43(3):525–546. doi: 10.1137/S0036144500378302. [Cross Ref]
  • Kandel ER, Schwartz J, Jessell T. Principles of neural science. 4. New York: McGraw-Hill; 2000.
  • Kramer MA, Kirsch HE, Szeri AJ. Pathological pattern formation and cortical propagation of epileptic seizures. Journal of the Royal Society Interface. 2005;2:113–127. doi: 10.1098/rsif.2004.0028. [PMC free article] [PubMed] [Cross Ref]
  • Kramer MA, Lopour BA, Kirsch HE, Szeri AJ. Bifurcation control of a seizing human cortex. Physical Review E. 2006;73:41928. doi: 10.1103/PhysRevE.73.041928. [PubMed] [Cross Ref]
  • Kramer MA, Szeri AJ, Sleigh JW, Kirsch HE. Mechanisms of seizure propagation in a cortical model. Journal of Computational Neuroscience. 2007;22(1):63–80. doi: 10.1007/s10827-006-9508-5. [PubMed] [Cross Ref]
  • Liley DTJ, Wright JJ. Intracortical connectivity of pyramidal and stellate cells: Estimates of synaptic densities and coupling symmetry. Network: Computation in Neural Systems. 1994;5:175–189. doi: 10.1088/0954-898X/5/2/004. [Cross Ref]
  • Liley DTJ, Cadusch PJ, Wright JJ. A continuum theory of electro-cortical activity. Neurocomputing. 1999;26–27:795–800. doi: 10.1016/S0925-2312(98)00149-0. [Cross Ref]
  • Liley DTJ, Cadusch PJ, Dafilis MP. A spatially continuous mean field theory of electrocortical activity. Network: Computation in Neural Systems. 2002;13:67–113. [PubMed]
  • Lopour BA, Szeri AJ. Advances in cognitive neurodynamics: Proceedings of the international conference on cognitive neurodynamics. Berlin: Springer; 2008. Chapter 86: Spatial considerations of feedback control for the suppression of epileptic seizures; pp. 495–500.
  • Nieuwenhuys R. The neocortex: An overview of its evolutionary development, structural organization and synaptology. Anatomy and Embryology. 1994;190(4):307–337. doi: 10.1007/BF00187291. [PubMed] [Cross Ref]
  • Nunez PL, Srinivasan R. Electric fields of the brain: The neurophysics of EEG. Oxford: Oxford University Press; 2006.
  • Richardson KA, Schiff SJ, Gluckman BJ. Control of traveling waves in the mammalian cortex. Physical Review Letters. 2005;94:28103. doi: 10.1103/PhysRevLett.94.028103. [PMC free article] [PubMed] [Cross Ref]
  • Robblee LS, Rose TL. Neural Prostheses: Fundamental Studies, chapter 2: Electrochemical guidelines for selection of protocols and electrode materials for neural stimulation. Englewood Cliffs: Prentice Hall; 1990.
  • Rubinstein JT, Spelman FA, Soma M, Suesserman MF. Current density profiles of surface mounted and recessed electrodes for neural prostheses. IEEE Transactions on Biomedical Engineering, BME. 1987;34(11):864–875. doi: 10.1109/TBME.1987.326007. [PubMed] [Cross Ref]
  • Schiff SJ, Sauer T. Kalman filter control of a model of spatiotemporal cortical dynamics. Journal of Neural Engineering. 2008;5:1–8. doi: 10.1088/1741-2560/5/1/001. [PMC free article] [PubMed] [Cross Ref]
  • Schiff SJ, Huang X, Wu J-Y. Dynamical evolution of spatiotemporal patterns in mammalian middle cortex. Physical Review Letters. 2007;98:178102. doi: 10.1103/PhysRevLett.98.178102. [PMC free article] [PubMed] [Cross Ref]
  • Spruston N. Pyramidal neurons: Dendritic structure and synaptic integration. Nature Reviews. Neuroscience. 2008;9(3):206–221. doi: 10.1038/nrn2286. [PubMed] [Cross Ref]
  • Steyn-Ross ML, Steyn-Ross DA, Sleigh JW, Liley DTJ. Theoretical electroencephalogram stationary spectrum for a white-noise-driven cortex: Evidence for a general anesthetic-induced phase transition. Physical Review E. 1999;60(6):7299–7311. doi: 10.1103/PhysRevE.60.7299. [PubMed] [Cross Ref]
  • Steyn-Ross ML, Steyn-Ross DA, Sleigh JW, Whiting DR. Theoretical predictions for spatial covariance of the electroencephalographic signal during the anesthetic-induced phase transition: Increased correlation length and emergence of spatial self-organization. Physical Review E. 2003;68:21902. doi: 10.1103/PhysRevE.68.021902. [PubMed] [Cross Ref]
  • Steyn-Ross ML, Steyn-Ross DA, Sleigh JW. Modelling general anaesthesia as a first-order phase transition in the cortex. Progress in Biophysics & Molecular Biology. 2004;85:369–385. doi: 10.1016/j.pbiomolbio.2004.02.001. [PubMed] [Cross Ref]
  • Suesserman MF, Spelman FA, Rubinstein JT. In vitro measurement and characterization of current density profiles produced by nonrecessed, simple recessed, and radially varying recessed stimulating electrodes. IEEE Transactions on Biomedical Engineering. 1991;38(5):401–408. doi: 10.1109/10.81558. [PubMed] [Cross Ref]
  • Sunderam, S., Chernyy, N., Mason, J., Peixoto, N., Weinstein, S. L., Schiff, S. J., et al. (2006). Seizure modulation with applied electric fields in chronically implanted animals. In Proceedings of the 28th IEEE EMBS Annual International Conference. [PMC free article] [PubMed]
  • The Epilepsy Foundation (2009). Treatment options: Surgery. World Wide Web,
  • Wilson MT, Steyn-Ross A, Sleigh JW, Steyn-Ross ML, Wilcocks LC, Gillies IP. The k-complex and slow oscillation in terms of a mean-field cortical model. Journal of Computational Neuroscience. 2006;21:243–257. doi: 10.1007/s10827-006-7948-6. [PubMed] [Cross Ref]

Articles from Springer Open Choice are provided here courtesy of Springer