## Abstract

Here we show that brain seizures can be effectively suppressed through random modulation of the brain medium. We use an established mesoscale cortical model in the form of a system of coupled stochastic partial differential equations. We show that by temporal and spatial randomization of parameters governing the firing rates of the excitatory and inhibitory neuron populations, seizure waves can be significantly suppressed. We find that the attenuation is the most effective when applied to the mean threshold potential. The proposed technique can serve as a non-invasive paradigm to mitigate epileptic seizures without knowing the location of the epileptic foci.

## 1. Introduction

Epilepsy, characterized by epileptic seizures, is the fourth most common neurological disorder in the world affecting over 50 million people [1]. Epileptic seizures start due to simultaneous electric-signal firing by a large number of nearby neurons. Details of the mechanism that causes this simultaneous firing are not yet fully understood, but when a bundle of neurons fire together, they collectively form a localized high-energy wave on the surface of the cortex. This localized high-energy wave then propagates over the surface of the brain, similar to how water waves propagate in a pond, and may result in a number of physical effects such as erratic shaking of the body, loss of consciousness, lethargy and pain.

Seizures occur unpredictably. Patients who suffer from epilepsy may be limited in their daily activity. The most common treatment, useful for approximately 70% of patients, is based on antiepileptic medications and special diets. The other 30% of patients may undergo invasive surgeries, such as removing seizing portions of the cortex or implanting neurostimulators in the body to alleviate their condition [2]. These invasive treatments do not guarantee complete suppression of epileptic seizures. Furthermore, the resection surgery risks damaging healthy brain tissue and can, therefore, cause loss of functionality. For these reasons, non-invasive treatments to attenuate epileptic seizures are highly desired.

Most existing methods for epileptic seizure control rely on targeting specific locations of the brain with electric or magnetic stimulation. These modalities often cannot be applied non-invasively. A number of open and closed loop feedback control methods have been proposed to mitigate epileptic seizures through transcranial electric or magnetic modulation [3,4]. Such methods require launching a lot of energy that may cause damage to the cortical tissue. Furthermore, the required resolution to target epileptic foci may not be easily achieved. Recently, optical stimulation has been proposed for modulation of epiletic activity based on optogenetic stimulation and inhibition of neural activity [5]. Attenuation of seizure activity through open-loop optogenetic control has also been computationally developed. By coupling a mesoscale optogenetic model with a mesoscale cortical model, a charge-balanced control scheme was introduced, minimizing damage to cortical tissue [2]. Optical stimulation of neurons requires expression of opsins in neurons through injection of virus vectors, a procedure which is not amenable for human testing.

*In vivo* experiments in rats have suggested that electric stimulation through feedback can attenuate seizure activity [6]. Attenuation of seizure activity by modulating the threshold potential with electric fields has been observed in rat brain slices *in vitro*. Closed and open loop feedback control of external electric fields for seizure modulation have been confirmed through computational models and simulations. Differential, linear and filter controllers were successfully implemented and were found to attenuate seizure activity [7]. A charge-balanced feedback control algorithm was implemented in the mesoscale cortical model to minimize cortical damage when applying electrical signals [8]. The use of control mechanisms serves as the underlying theme in all of these methods to modulate seizure activity. These non-invasive methods require sensors to determine when a seizure is occurring, process the data through a feedback system and then apply the proper signal to attenuate the seizure. However, proper stimulation with the correct feedback controller is crucial. This is important because feedback may not only fail to mitigate the seizure, but it can also amplify the seizure, disrupt normal brain activity and be detrimental to the patient. For example, it is shown that positive proportional feedback can amplify seizure activity, thus making the condition worse [6].

Here, we present another approach to seizure mitigation inspired from a wave phenomenon called Anderson localization. Localization is the dampening of waves by random modulation of the medium. For example, ocean waves lose their energy as they propagate over irregular topography (random medium). The concept was first discovered by P. W. Anderson in the context of solid-state physics, for which he won a Nobel Prize in 1977 [9]. Localization has since been observed in many other wave systems, and is used today for engineering purposes [10,11]. For example, in tsunami mitigation, one approach is to plant coastal vegetation (e.g. random trees) on the shoreline.

Epileptic seizures are due to collective pathological activity of a large population of neurons. Therefore, from a mesoscale perspective, a seizure can be considered as a synchronized collective of waves propagating across cortical regions. We show here that this wave propagation behaviour can be disrupted through local modulation of neurons (hypothetically, through a phased-array of ultrasound transducers) to create multiple local modes that can dissipate the seizure wave energy. A major advantage of the proposed methodology is that we do not need to target the epileptic foci. Therefore, the detection mechanism does not require high resolution or a dense array of electroencephalogram (EEG) sensors. We only need to monitor the magnitude of the epileptic seizure wave (for example, using EEG) and disrupt it by random modulation of neural activity in the cortex.

While details of how the required modulation is physically implemented are beyond the scope of this paper and will be reported elsewhere, we will discuss the possible use of ultrasonic waves to achieve this purpose. It has been recently demonstrated that ultrasonic waves can modulate (stimulate and inhibit) neural activity [4]. Ultrasound is a medically safe imaging modality and has been used for imaging internal organs and fetuses. In the treatment of brain cancer, high-frequency ultrasound is used to monitor apoptosis in malignant tumours [12]. It is conceivable that ultrasound can be used to trigger mechano-signal transduction pathways in biological tissues and cells. Ultrasound in certain frequencies (0.6–0.7 MHz) has minimal propagation loss in the brain (less than 0.5 dB cm^{−1} MHz), has maximal transcranial transmission, does not interfere with electromagnetic waves and has a relatively small wavelength [4]. First use of ultrasound for modulating neural activity was reported by Fry *et al.* [13–16].

In this paper, we discuss the concept of Anderson localization for suppression of tonic–clonic epileptic seizures in the context of a cortical model. Our proposed method may be implemented in practice using a variety of approaches including electrical, magnetic, ultrasonic or optical stimulation. Although there are different types of seizure activities such as absence, myoclonic or atonic seizures, we focus on the generalized tonic–clonic seizure activity. The cortical model we use has already been studied and verified to model this type of seizure [1,2].

## 2. Mesoscale cortical model

The cortex is a complex network of interconnected neural circuits consisting of billions of neurons. Precise modelling of all neural interactions in the human cortex is an elusive goal of neuroscience and is computationally intractable. In the past two decades, a tractable mesoscale model of the cortex in the form of a set of stochastic partial differential equations (SPDEs) has been introduced [17,18]. Rather than attempting to model every neuron and neural connection in the cortex, the equations model the mesoscale electrical dynamics in the cortex. That is, instead of examining the dynamics of the soma membrane potential of each individual neuron, the mesoscale model considers the spatial average of the potential of a population of neurons. The mesoscale model is expressed as a set of equations that includes anatomically derived neuronal connectivities, ionic reversal potentials, and fast excitatory and inhibitory channel kinetics.

The excitatory mean soma potential *h*_{e} that is modelled by the mesoscale equations are related to electrocorticogram (ECoG) and EEG recordings. Electrodes used in EEG and ECoG record electrical activity of a local region of the cortex. The data obtained from these recordings are what the mesoscale equations aim to model. In fact, both EEG and ECoG are related to the spatially averaged local field potential (LFP), which is measurable on the scalp (in the case of EEG) or on the surface of the cortex (in the case of ECoG). In [18], the author argues that the deviation between the excitatory mean soma potential and the resting potential is proportional to the LFP. As EEG data are linearly related to the spatially averaged LFP, we can conclude that EEG data are proportional to the spatially averaged excitatory soma potential.

This mesoscale model has been successfully used to study phenomena such as epileptic seizures [1,2,7,8,19], sleep [20,21] and anaesthesia [22–24]. In [1], the authors compared time series data of real seizures with oscillatory activity exhibited by the mesoscale cortical model. More specifically, they compared the peak frequency and wave speed of the observed seizure data with computational simulations from the model. Therefore, Kramer *et al*. [1] validates our use of this model to computationally produce seizures.

The dimensionless mesoscale cortical model relates (i) the excitatory and inhibitory soma potentials and _{i}, (ii) the postsynaptic activation input from local excitatory and inhibitory neurons to excitatory neurons *Ĩ*_{ee} and *Ĩ*_{ie}, (iii) the postsynaptic activation input from local excitatory and inhibitory neurons to inhibitory neurons *Ĩ*_{ei} and *Ĩ*_{ii}, and (iv) the long range inputs from excitatory neurons to excitatory and inhibitory neurons and , through the following coupled SPDEs:
2.1*a*
2.1*b*
2.1*c*
2.1*d*
2.1*e*
2.1*f*
2.1*g*
2.1*h*The subscripts ‘e’ and ‘i’ refer to the excitatory and inhibitory neuron populations, respectively. The double subscripts such as ‘ei’, for example, mean from the ‘e’ to the ‘i’ population. For instance, *Ĩ*_{ie} represents postsynaptic activation input from inhibitory neurons to local excitatory neurons. The model consolidates neurons into two types, spatially averaged excitatory and inhibitory neurons, in which local interactions with glial and extracellular matrix are also incorporated. That is, neurons and the media surrounding the neurons are merged into a single entity called the spatially averaged neuron. The connections between the neurons are consolidated into local (neighbouring) and long range (cortico-cortical) ^{1} interactions.

Equation (2.1*a*) describes the time evolution of the excitatory soma potential _{e} as a function of itself and the postsynaptic activation input from local excitatory *Ĩ*_{ee} and inhibitory neurons *Ĩ*_{ie}. Similarly, equation (2.1*b*) describes the time evolution of the inhibitory soma potential _{i} as a function of itself and the postsynaptic activation input from local excitatory *Ĩ*_{ei} and inhibitory neurons *Ĩ*_{ii}. Equations (2.1*c*) through (2.1*f*) relate the time evolution of the postsynaptic activation inputs between local neurons as functions of firing rates of neurons , long range cortico-cortical inputs from excitatory neurons , and subcortical inputs ,, where are the stochastic components.

Equations (2.1*g*) and (2.1*h*) relate the time and space evolution of the long range cortico-cortical inputs with the firing rate of the excitatory population. Each of the eight dynamical variables are functions of dimensionless space and time (). The topology is implicitly defined by the equations. A graph of the relationships between all these connections can be easily drawn from the equations. For a diagram of all these connections, see [18].

The two firing rate functions
2.2*a*and
2.2*b*model the firing rates of the exictatory and inhibitory neuron populations, respectively, as a function of their membrane potentials. These sigmoid functions were originally derived in [18] and then reformulated to their dimensionless form in [1]. We discuss these functions in further detail in §3.2. In this paper, we will modify these functions to incorporate coupling between the parameters in the functions. A table defining all the coefficients in this model is provided in appendix A.

## 3. Implementation

### 3.1. Numerical details

We now describe how the mesoscale cortical model is numerically implemented for simulations. The numerical method of lines [25] is adopted to solve the system of SPDEs. The spatial derivatives are discretized first to create a system of ordinary differential equations which is integrated forward in time. The second derivative in space is discretized with a second-order accurate central difference formula and the two-dimensional Laplacian operator is discretized using a 5-point stencil [25]. The resulting system of ordinary differential equations (with stochastic forcing) are integrated forward in time with a fourth order Runge–Kutta method.

We now explain how the stochastic inputs in the equations are implemented. Suppose the spatial step size is Δ, and the time step size is . Random inputs at location and time are determined by drawing from the following independent normal random variables [8]:
3.1*a*
3.1*b*
3.1*c*
3.1*d*where are independent, identically distributed standard normal random variables, and are user-determined parameters. We discuss how the value *α* is chosen in §3.1.1.

To explain what is happening, let us only consider . The same reasoning will apply to . Recall from (2.1) that the subcortical input from excitatory neurons to excitatory neurons is , where *P*_{ee} is a *constant* and is a *stochastic* input. This modelling choice is saying that the subcortical input is a random process in which at every point in space and time, the value of the subcortical input is determined randomly by drawing from some random variable. A Gaussian white noise process is a common choice to model [1,2,8] because when it is integrated in time, the result is a Brownian motion, and because the fluctuations are uncorrelated (*δ*-function correlated) [26]. As we discretize the cortex to simulate the model, we must also discretize the white noise process. We do this by drawing from independent, identically distributed standard normal random variables for each space and time discretization point (*m,n*), and then scaling the variance by multiplying each sample by . The variance is scaled according to the spatial and temporal discretizations so that the noise is *δ*-function correlated as . Note that one can control the variance of the subcortical inputs by changing α_{ee}.

#### 3.1.1. Initial and boundary conditions

The initial conditions we choose are related to the equilibrium solution of the system, which is the solution to the model (without stochastic inputs) that does not change in space or time. The equilibrium solution is found by setting the time and spatial derivatives to be zero, and the stochastic inputs to be zero (, , ). This will result in a set of algebraic equations for which we can solve for the state variables. Setting the initial conditions precisely at equilibrium is unrealistic as the cortex is never in a uniform equilibrium state. Instead, we add noise to the equilibrium conditions to get a more realistic initial condition. So for every point in space, we offset the equilibrium values by randomly adding a value drawn from a normal distribution. The addition of random noise is an attempt to model normal brain activity before it goes into a seizure state.

Whether or not a seizure forms in this model is mostly dependent on the influence of input on the excitatory soma potential (*Γ*_{e}) and subcortical input to the excitatory neuron populations (*P*_{ee}). This hypothesis is thoroughly argued in [1]. The parameter *α*, which controls the magnitude of randomness from subcortical inputs, also matters. However, there is neither a commonly accepted value nor a physically meaningful way for determining *α*. For example, Selvaraj *et al.* [2] assumes a large range of 0.001–1.15 for *α* to demonstrate the robustness of optogenetic control. In the upcoming one-dimensional simulations of the model, we chose *α* so that there would be some normal brain activity in the first second before exhibiting seizure activity (which is manifested in the form of oscillations in space and time). This can be seen in figure 1. If the initial conditions are too close to equilibrium or if the noise in the system *α* is too weak, then seizure activity may not begin for a long time.

Periodic boundary conditions are implemented for both the one- and two-dimensional versions of the cortical model. This is a common modelling assumption as we do not know the exact geometry of the cortex [1,2].

One may ask how does one know if the data produced from the chosen initial conditions correspond with seizures. Recall from §2 that in [1], the authors compared time series data of real seizures with simulated data generated from the mesoscale model. They found that the average frequency of maximum power and wave speed of the simulated seizure matched well with real seizures. Therefore, we may safely assume that the results from this model accurately simulate seizures.

#### 3.1.2. Examples of simulated seizures

We show examples of seizure formation and propagation in figures 1 and 2. The main quantity we are studying is the excitatory soma potential _{e}. In the one-dimensional case (figure 1), the spatial step size is , and the temporal step size is Δ = 0.0025. These correspond to dimensional step sizes of Δ*x* = 14 mm and Δ*t* = 0.1 ms. As neurological processes take place with the resolution of about 1 ms (approximately equal to the refractory period), it is valid to use time steps of around that magnitude to properly capture the dynamics of the system. Furthermore, using smaller spatial time steps may be invalid as the model we are using is a mesoscale model. More specifically, as the mesoscale model is related to the LFP, there would be no valid physical justification to use an arbitrarily small spatial step. Most parameters chosen for the one-dimensional case are the baseline cortical parameters described in [1], which are listed in table 2 in appendix A, except for . These amended parameters correspond to a seizure state of the cortex (cf. [1]).

In figure 2, we show a two-dimensional seizure simulation. We use the same parameters as the one-dimensional case except for *α* = 0.002 and the temporal step size of . This step size corresponds to the dimensional step size of Δ*t* = 1 ms (cf. [2]). Figure 2 shows that even though the initial conditions are chosen such that the cortex is near equilibrium, the stochastic inputs lead to strong coherent oscillations which correspond to seizures. The more interesting result is the propagation of spiral waves through space. This further demonstrates and supports the wave like nature of seizures as presented by past computational research [2,27]. While the actual surface area of the human cortex is around 500 × 500 mm, we can better demonstrate the wave nature developed by seizure waves by adopting the larger 1400 × 1400 mm domain [2].

### 3.2. Modifying the firing rate functions

Now we come back to discuss the neuron firing rate functions, which will become important for localization. We first restate equation (2.2) in its *dimensional* form to get a clearer picture of all the parameters involved in this model:
3.2This model is sufficient in describing the firing rate of each population of neurons as a function of their soma potentials when the parameters governing the firing rate, *g*_{j}, which describes the spread of threshold potentials across neurons, *θ*_{j}, which is the population mean threshold potential and *S*^{max}_{j}, the maximum firing rate, are held constant. However, this equation is not sufficient when the parameters of the function may be changing in space or time. We may naively assume that to incoporate spatial and temporal inhomogeneity in the parameters, we only need to replace the constants with functions. This will not work because this model is only an approximation of a more complex model that may no longer be valid with inhomogeneous parameters.

Equation (3.2) is an approximation of the more accurate firing rate function of type neuron population
3.3where *r*_{s} is the absolute refractory period, the time in which a neuron is firing and obviously cannot be triggered to fire again. It is common to assume that such that equation (3.3) reduces to equation (3.2) [18]; however, as we are considering randomizing the parameters (i.e. letting them vary with time and space), neglecting this value can lead to wrong results. Equation (3.4) shows how the maximum firing rate is coupled with the threshold potential. This equation can be derived by finding the inflection point of (3.3). In the upcoming simulations, we assume *r*_{s} = 1 ms, which is about the time for a neuron to be depolarized past the activation potential threshold, fire, and then polarize [18]. From appendix A, *S*^{max} = 500 Hz so , which is not negligible. Equation (3.4) shows coupling between *S*_{j}^{max}, *θ*_{j} and *g*_{j}. Fixing any one of the three parameters, while changing another parameter will invariably change the third parameter. That is, no one parameter can be changed independently of the other two parameters. Clearly we see that randomizing parameters in the approximation (3.2) will not exhibit this coupling.
3.4With the incorporation of the *r*_{s}*S*_{j}^{max} term, we need to calculate new *θ*_{e} and *θ*_{i} values as they were offset by and , respectively. The new dimensionless firing rate function is therefore obtained as
3.5Note that we do not change as it is already a dimensionless term. All the upcoming simulations use the modified dimensionless firing rate function.

### 3.3. Implementing localization

To disrupt seizure waves, we apply the Anderson localization idea here. This requires randomizing the cortical medium through which the seizure wave is propagating. This is a general concept and can be applied to any type of wave system. For example, in the case of ocean waves, the seabed is randomized by small ripples to achieve localization [10,11]. For the electrical waves propagating on the cortex, table 2 in appendix A presents all the possible parameters in the scope of the mesoscopic model that can be potentially tweaked to disrupt the seizure. Most of these parameters are difficult to change externally because they characterize the anatomical and physiological structure of the cortex. For example, changing parameters such as the number of synaptic connection between neurons , or the sub-cortical inputs *P*_{ee,ei,ie,ii} to the neurons is formidable as they define physiological and deep intercortical interactions and cannot be changed easily. We, therefore, focus our attention on the parameters governing the firing rates of the populations of neurons. The physical mechanism of randomizing these parameters is beyond the scope of this paper. However, a few possible methods include using non-invasive ultrasound, electrical stimulation or optogenetic stimulation to modulate the firing rates of neurons [4].

To achieve localization, the constant parameters *S*^{max}_{j}, *g*_{j} or *θ*_{j} will be replaced with random processes, or random functions of space and time. How this is implemented is as follows:
3.6
3.7
3.8where *σ*_{k} are the standard deviations of the normal distribution . In other words, at each point in space and time, we replace the original parameter, say *θ*_{e}, with a sample drawn from a normal distribution with mean *θ*_{e} and variance . As there are two populations of neurons, there are a total of six parameters that can be randomized: *S*^{max}_{e}, *S*^{max}_{i}, *g*_{e}, *g*_{i}, *θ*_{e} and *θ*_{i}. In the following section, each of these parameters will be randomized over a range of variances to see how they affect the propagation of the seizure waves. We begin with a discussion of simulations that will serve as the control cases. Then the effectiveness of localization is first discussed for the one-dimensional model, which will motivate discussion and provide insight for the two-dimensional model.

## 4. Seizure suppression in one-dimensional mesoscale model

The one-dimensional cortical model with periodic boundary conditions can be interpreted as modelling EEG or ECoG data with linear strip electrodes around the scalp [1]. In this section, we will simulate seizures using the one-dimensional cortical model and suppress them via localization. Randomizations in both space and time are run for each of the six parameters over a range of variances. We first allow the seizure to form in the first three seconds of the simulations. Afterwards, localization is applied for a period of three seconds. The vertical line at *t* = 3 s in the following plots shows the onset of applying the localization. The control case (i.e. with no localization applied) is shown in figure 1.

Each simulation uses the exact same stochastic inputs and initial conditions. The spatial step size is *h* = 0.05 and the temporal step size is *k* = 0.0025. These correspond to dimensional step sizes of Δ*x* = 14 mm and Δ*t* = 0.1 ms. All results are confirmed to be converged. Specifically, we tested a few cases of localization with different spatial and temporal step sizes and showed that the effects of localization are not dependent on the step size. In the discussion, we provide a brief demonstration of convergence.

### 4.1. Qualitative effects of localization

We begin by showing simulations of localization in action, and discuss their qualitative features.

#### 4.1.1. Randomizing *θ*_{e}

Randomizing *θ*_{e}, the mean excitatory threshold potential in space and time suppresses seizures successfully by reducing the seizure amplitude. Figure 3*a*,*b* shows localization applied through randomizing *θ*_{e} in both space and time with *σ* = 0.20. We can see that the seizure is completely attenuated within 0.5 s throughout the domain due to the random temporal and spatial perturbations. Note that the oscillations are centred around the equilibrium potential *h*_{e}∼−53 mV. This result is similar to the closed loop control of an external electric field to suppress seizure activity [8].

#### 4.1.2. Randomizing *θ*_{i}

Randomizing *θ*_{i}, the inhibitory population threshold potential decreases the frequency of the seizure oscillations as the underlying randomization increases in variance. Figure 4 shows localization being applied by randomizing *θ*_{i} in both space and time. The general trend is that low values of *σ*, such as *σ* = 0.15 (figure 4*a*,*b*), show little to no seizure attenuation at all. However, as *σ* is increased, we see seizure attenuating in an unusual way. Rather than suppressing the amplitude of the oscillations in space and time, the frequency of the seizures is reduced. In figure 4*c*,*d*, *σ* = 0.30, we see that the seizure is more spaced out, suggesting that the frequency has decreased. This is confirmed in figure 4*d* in which we can see that the frequency of the seizure has decreased dramatically, but the amplitude is only slightly affected. If *σ* is increased further, we will see that the frequency will be further reduced until the seizure is completely eliminated (not shown here). Observe that oscillations in excitatory soma potential are centred around the resting potential of the neuron, *h*_{e}∼−70 mV. This result is similar to the open loop optogenetic control of inhibitory neurons [2].

#### 4.1.3. Randomizing *S*_{e}^{max}

Randomizing *S*^{max}_{e}, the maximum excitatory firing rate, will also decrease the amplitude of seizures, although only with higher variances than what was needed for *θ*_{e}. Figure 5*a*,*b* shows localization being applied by randomizing *S*^{max}_{e} in both space and time with *σ* = 0.3. The mechanism for the suppression is most probably due to the fact that it is coupled with *θ*_{e}. That is, randomizing *S*^{max}_{e} also indirectly randomizes *θ*_{e}. But we see some sort of ‘delay’ as the seizure does not attenuate as quickly or as uniformly. In the *θ*_{e} cases, we saw that seizure attenuation started soon after *t* = 3 s and that the seizure attenuated uniformly everywhere.

#### 4.1.4. Randomizing *S*_{i}^{max}

Randomizing *S*_{i}^{max}, the maximum firing rate of the inhibitory population does not sufficiently attenuate seizures. Figure 5*c*,*d* shows localization being applied by randomizing *S*^{max}_{i} with *σ* = 0.50 in both space and time. While extremely subtle, the frequency of the seizures is, in fact, being decreased as *σ* is increased. The frequency before localization is approximately *f*_{0} = 9.5 Hz, and after is *f*_{1} = 7 Hz. As we are applying localization to inhibitory neurons, this effect is consistent with the case when applying randomization to *θ*_{i}, but with less success. Even with a larger amplitude of noise, i.e. *σ* = 0.50, we can barely see any attenuation.

#### 4.1.5. Randomizing *g*_{e}

Randomizing *g*_{e} will attenuate seizure amplitude, but only with high variances. Figure 6 shows localization being applied by randomizing *g*_{e} in space and time with *σ* = 0.25 and *σ* = 0.5. This example shows that randomizing *g*_{e} with low variances does not have a noticeable effect on the seizure; only at high variances do we see seizure attenuation happen relatively quickly throughout the domain. It is important to note that the level the standard deviations needs to be to attenuate seizure for *g*_{e} is much higher than for other parameters.

#### 4.1.6. Randomizing *g*_{i}

Randomizing *g*_{i} reduces seizure frequency, but does not eliminate seizures completely. Figure 7 shows localization being applied through randomizing *g*_{i} with *σ* = 0.25 shown in figure 7*a*,*b*, and *σ* = 0.5 shown in figure 7*c*,*d*. This figure demonstrates the difficulty of attenuating seizures by randomizing *g*_{i}. As we are applying localization to inhibitory neurons, we see that the frequency of the seizures is decreased as *σ* is increased, which is similar to the cases for *θ*_{i}. However, it requires much higher *σ* to achieve similar levels of attenuation.

### 4.2. Evaluating efficacy of localization quantitatively

From figures 3–7, one can roughly see that by randomizing and *S*_{e}, the oscillatory behaviour of a seizing cortex is suppressed. To quantify and compare the efficacy of localization, we need to define a metric for seizure suppression.

Typically for wave phenomena, one looks at the frequency and amplitude of the wave to estimate the energy. In our case, looking at these quantities alone will not be a sufficient measure because we have seen that the seizures are attenuated in two ways: a decrease in amplitude or a decrease in frequency. Therefore, a better metric must be used. In [28], the authors propose a metric called *ictality*, which is based on the autocovariance of the excitatory soma potential over time at a single point in space. Given data of the excitatory soma potential _{e} over time at some point in space, we can calculate the ictality by first plotting the autocovariance of _{e}, and then dividing the height of the second-order peak by the first-order peak. This procedure is equivalent to plotting the autocorrelation of _{e} and finding the value of the second peak as the autocorrelation is the autocovariance normalized according to the variance of the signal. See figure 8 for an example of calculating the ictality of a seizure.

Ictality measures how periodic _{e} is over time. If the signal is close to a sinusoid, then the data is highly correlated with itself, so the autocorrelation will decay very slowly, and the ictality will be high. If the signal is very noisy with no distinct period, then the data will not be correlated, so the autocorrelation will decay very quickly, and the ictality will be low. Seizures are characterized by oscillatory behaviour of _{e}. A cortex undergoing a seizure will see _{e} be more correlated with itself. Therefore, during a seizure, the ictality will be high. We can use this measure to see if ictality of the seizure decays with greater randomization of the parameters.

As localization relies on random inputs, a more statistical study is required. We perform an ensemble of simulations for the one-dimensional case. For each parameter ( ), we performed approximately 10 simulations for each relative standard deviation from *σ* = 0.05 to *σ* = 0.5. Then for each simulation, we calculated the ictality of the seizure 0.5 s *after* localization is applied at every point in the domain and took their mean. These data are plotted in figure 9.

Clearly this shows that the best candidates for localization are *θ*_{e} and *S*_{e}. While *θ*_{i} and *S*_{i} also show decay, these results are questionable as the attenuated seizure converges near the resting potential rather than at the equilibrium potential. Furthermore, this also shows that the parameters *g*_{e} and *g*_{i} are not effective for mitigating seizures, or at least far less effective than other parameters.

## 5. Seizure suppression in two dimensions

We now apply localization for all of the six aforementioned parameters by randomizing them at different variances in the two-dimensional model. We will compare the upcoming simulations to the control simulation in figure 10, which was run for total time *t* = 1 s on a 700 × 700 mm size domain. The parameters used to run these simulation are the baseline parameters of the seizure state as defined in [1] and used in [2] (except for *Γ*_{e} = 0.66 × 10^{−3}, *P*_{ee} = 548.00, *α* = 1.15) in order to compare effectiveness of seizure attenuation with past methods such as optogenetic control [2] and closed loop feedback through electrical stimulation [8]. Periodic boundary conditions are implemented and the initial conditions are the equilibrium solution that is randomly perturbed. For each of the following simulations, we provide two plots. Plot (*a*) shows a one-dimensional cross section of the plane over time at *y* = 350mm. Plot (*b*) shows the evolution of the excitatory soma potential at a single point in the one-dimensional cross section at *x* = 350 mm. Localization is applied at *t* = 0.5 s for all simulations. The spatial step size is *h* = 0.05, which corresponds to 14 mm for both x- and y-axes. The temporal step size is *k* = 0.025, which corresponds to 1 ms.

### 5.1. Simulations

#### 5.1.1. Randomizing *θ*_{e}

Similar to the one-dimensional case, randomization of *θ*_{e} reduces the amplitude of seizures. Figure 11*a*,*b* shows localization applied to *θ*_{e} with *σ* = 0.25 through randomization in both space and time. Randomizing the medium interrupts the seizure propagation and attenuates the seizure. Qualitatively, we note that the attenuation in two dimensions is not as effective as in the one-dimensional case. The residual oscillations in the two-dimensional case have larger amplitudes.

#### 5.1.2. Randomizing *θ*_{i}

Randomizing *θ*_{i} suppresses seizures by reducing the seizure wave frequency. In figure 12*a*,*b*, localization is applied by randomizing *θ*_{i} in both space and time. Note that it is far easier to attenuate seizure by modulating the threshold potential of inhibitory neurons in the two-dimensional case than in the one-dimensional case. In the two-dimensional case, randomizing *θ*_{i} with *σ* = 0.15 will result in complete attenuation while in the one-dimensional case, we need σ = 0.45. Similar to the one-dimensional case, random modulation of the inhibitory neurons decreases the frequency rather than the amplitude of the seizures.

#### 5.1.3. Randomizing *g*_{i}

Figure 12*c*,*d* shows localization being applied by randomizing *g*_{i} in both space and time with *σ* = 0.2. Note that as we are modulating inhibitory neurons, the excitatory soma potential converges near the resting potential of neurons, rather than near equilibrium. Fluctuating near the resting potential is not representative of a real life scenario. Also note that the seizure is completely eliminated, which contrasts with the one-dimensional simulations. Figure 7 shows that the seizures are beginning to attenuate, but at *σ* = 0.5, the seizure is barely perturbed. However, figure 12*c*,*d* shows that randomizing *g*_{i} with just *σ* = 0.2, we can see complete attenuation.

#### 5.1.4. Randomizing *g*_{e}, *S*_{e}^{max}, *S*_{i}^{max}

Randomizing *g*_{e} does not attenuate seizures. Figure 13*a*,*b* shows localization being applied by randomizing *g*_{e} in both space and time with *σ* = 0.3. This confirms the one-dimensional result that randomizing *g*_{e} to attenuate seizures is quite difficult. The amplitude decreases very slightly after localization. Also, from figure 13*c*–*f*, we see that randomizing *S*_{e}^{max}, *S*_{i}^{max} does not attenuate seizures in two dimensions. In neither case do we see much seizure attenuation as compared to the control case. This is different from the one-dimensional modulation of these parameters, where we saw some attenuation at *σ* = 0.25.

## 6. Discussion

Using a mesoscale model of the cortex, described by the set of coupled nonlinear SPDEs in §2, we showed how randomly modulating different parameters governing the firing rates of neurons affects the propagation of seizure waves in the brain. From the results presented in §§4 and 5, we can see that localization of seizure waves is the most effective when randomization is applied to *θ*_{e} and *θ*_{i}, the parameters describing the population mean threshold potential for excitatory and inhibitory, respectively. In the case of randomizing *θ*_{e}, using a standard deviation of only *σ* = 0.20 times *θ*_{e}, reduced the seizure to incoherent oscillations (figure 3). For *θ*_{i}, we saw annihilation of the seizure, in which the frequency was reduced (figure 4). We achieved similar results in two dimensions for *θ*_{e} when *σ* = 0.25 (figure 11) and for *θ*_{i} when *σ* = 0.15 (figure 12*a*,*b*). However, there is a key distinction between the two cases: randomization of *θ*_{e} results in _{e} converging close to the equilibrium state of the cortex, while randomization of *θ*_{i} results in _{e} converging to a fixed state of the system first described by [1] in dynamical systems theory, which happens to be approximately the resting potential of neurons. The latter case cannot be a real physiological state of the cortex after seizure disruption. Another interesting difference is that randomizing *θ*_{e} reduces the seizure amplitude, but not the frequency (figures 3 and 11), while randomizing *θ*_{i} reduces the seizure frequency, but not the amplitude (figures 4 and 12*a*,*b*). These phenomena can be explained by looking at how excitatory and inhibitory neurons function in cortical activity. Among all the possible variables that we could modulate in the mesoscale model of the cortex, we found these two variables to be the most effective in attenuating the electrical activity in the fastest possible way. It is quite conceivable that randomizing a proper combination of theses two variables might be the optimum scenario to localize the seizures, which of course needs to be tested against practical limitations.

In the one-dimensional model, randomizing *g*_{e} and *g*_{i} resulted in slight attenuation of seizures at the levels where suppression occurred for the *θ*_{j} cases. At *σ* = 0.25 (figures 6*a*,*b* and 7*a*,*b*), there is very little observable attenuation. Only at *σ* = 0.5 (figures 6*c*,*d* and 7*c*,*d*) do we see attenuation for the *g*_{e} case and some semblance of efficacy for the *g*_{i} case. In two dimensions, randomizing *g*_{e} results in little attenuation (figure 13*a*,*b*). For *g*_{i}, we find that with *σ* = 0.2, seizures are also completely eliminated (figure 12*c*,*d*). However, this result is questionable as it was not as prominent in the one-dimensional case. The efficacy of randomization of the maximum firing rate parameters is harder to judge. In one dimension, we saw that randomization of *S*_{e}^{max} attenuated seizures for *σ* = 0.3 (figure 5*a*,*b*), but no noticeable effect could be seen for *S*_{i}^{max} at *σ* = 0.5 (figure 5*c*,*d*). In two dimensions, randomizing neither *S*_{e}^{max} nor *S*_{i}^{max} attenuated seizures (figure 13*c*–*f*).

To quantify how effectively we can mitigate seizure waves, we need to define a figure of merit. Defining a universal figure of merit is challenging as modulating the excitatory threshold potential will reduce seizure amplitude, while modulating inhibitory neuron threshold potential will reduce seizure wave frequency. Therefore, a measure simply based on the amplitude will fail to capture the mitigation of seizure in the latter case. The authors in [28] argue that looking at the autocorrelation plots of _{e} provides a good metric for intensity of seizures, in which the magnitude of seizure strength is related to how fast the autocorrelation plot of a seizure at a point decays. Therefore, their metric was applied to large batches of simulations and was found to be a good measure of seizure strength. We are able to quantitatively see in figure 9 that the ictality of the seizure decays with *σ* and that *θ*_{e} is the best candidate to apply localization.

Recall that with a modified firing rate function, we established that these model parameters are coupled with each other. Therefore, randomization of one parameter would result in the randomization of other parameters. Such couplings are usually ignored. However, we established coupling between the maximum firing rate and the population threshold potential and showed that completely neglecting the coupling will be insufficient in applying localization to the model. For future work, we can improve the model by considering deterministic spatial and temporal inhomogeneities of these parameters, which can cause further possible coupling effects. This would require re-derivation of the population neuron firing rate transfer function.

In all of the results we presented, the cortex was approximated as discrete points and randomization was applied at each grid point. We show here in a few exemplary cases that our results are not grid dependent and are converged. Figure 14 shows two cases of randomization in both space and time being applied to *θ*_{e} with *σ* = 0.2 in the two-dimensional cortical model. In figure 14*a*, we use Δ*x* = 14 mm and Δ*t* = 0.1 ms. In figure 14*b*, we use Δ*x* = 7 mm and Δ*t* = 7 ms. Note that the underlying randomness is not the same. These show similar effectiveness of randomization. The former case is a reduction of the time step by a factor of 10, the latter case is a reduction of the spatial step by two. These plots show that the effectiveness of localization is not dependent on grid size or underlying randomness and are indeed converged.

Further research is needed to translate the randomization of mesoscale model parameters to physical means of modulating the neural activity in the cortex. For example, it has been shown that cortical neural circuits can be stimulated using non-invasive ultrasound waves [4], where both excitatory and inhibitory pathways can be stimulated by using the ultrasound at the appropriate frequency. Also, it has been shown that the neural firing rate in the retina can be modulated by using ultrasound at different frequencies or different powers [29]. The implementation of the intervention mechanism can be based on electrical, magnetic or ultrasonic stimulation of local neural activity. We propose that an array of ultrasonic transducers can be designed to generate the patterns of interest in the cortex to localize seizure waves. It should be noted that in this method, we would not need prior knowledge of the epileptic foci and the ultrasonic array would generate the ‘random’ pattern of interest to disrupt epilepsy seizure waves. This would become even more interesting when the ultrasonic array would function in an adaptive way to account for variations among subjects. Obviously, such an implementation requires careful study of potential side effects on cortical regions.

In the proposed technique, seizures would be disrupted when the randomization system is on. Therefore, we are not claiming that the parameters of the biological system would be permanently changed. The spontaneous phase transitions as a result of an external intervention can prevent seizure waves from propagation throughout the brain right after the onset of seizure. As a future clinical utility of this technique, one could imagine a wearable device that would launch a pattern of ultrasound waves from outside the skull to disrupt seizure. Having said that, we cannot rule out the potential of long-term permanent changes that might result in curing epilepsy through a gradual *in situ* learning process.

Finally, we comment that a seizure state has been shown to be correlated to a Hopf bifurcation in the phase space of *h*_{e}−*Γ*_{e} parameters of the mesoscale model [1]. In other words, the brain goes into a seizure state when a stable equilibrium in *h*_{e}−*Γ*_{e} phase diagram turns into an unstable equilibrium and two new stable equilibria are formed about this newly formed unstable equilibrium. Our simulations show that with the introduction of proper excitations, the Hopf bifurcation in the phase diagram no longer exists.

## 7. Conclusion

In this paper, we have shown the possibility of using Anderson localization to disrupt seizure waves exhibited by a mesoscale cortical model. Borrowing the concept of Anderson localization from solid state physics, we have presented here, for the first time, the effects of random modulation on the cortical model as a means of seizure attenuation by randomizing the parameters governing the firing rate of a population of neurons. We first showed the control cases for comparison and confirmation of past research. The one-dimensional cortical model characteristics described in [1] and two-dimensional cortical model characteristics [2] were considered. The travelling waves in the one-dimensional model and the spiral like waves described by [2,27] were confirmed. We also took into account possible coupling between parameters, and used a more precise equation to describe the firing rates of the two populations as a function of their soma potentials in our simulations to properly explore localization. We demonstrated the possibility of seizure suppression by randomizing the six parameters affecting the firing rate transfer function in both one- and two-dimensional cases. Spatial and temporal randomization were applied to each of these parameters at different variances to observe the effect over a large range of randomization parameters. We identified the most effective parameters to be the population mean threshold potential for excitatory and inhibitory pathways, i.e. *θ*_{e} and *θ*_{i}. In our proposed method, disruption of seizure waves can be simply carried out using random modulation and does not require any knowledge of how exactly seizure waves are formed in the cortex. We hope that this computational study will motivate future experimental and mathematical studies.

## Competing interests

We declare we have no competing interests.

## Funding

The authors are grateful for the kind support from the Hellman Foundation.

## Appendix A

Tables 1 and 2 define and summarize the variables and parameters of the cortical model. The dimensional excitatory soma potential *h*_{e} is related to its dimensionless counterpart _{e} by , where is the resting potential. [1]. Similarly, dimensional space and time are related to their dimensionless form by and , where is the membrane time constant and is the mean axonal conduction speed. It is now established that the formation and propagation of seizures are due to abnormally high levels of subcortical input to the excitatory neuron population *P*_{ee} and *Γ*_{e} [1,19].

## Footnotes

↵1 Technically should be written as , respectively. But as the model does not consider , the first subscripts are dropped.

- Received October 31, 2016.
- Accepted January 16, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.