Brain model state space reconstruction using an LSTM neural network

Scalp and intracranial EEG record electrical activity signals captured on the scalp and on the surface of the brain, respectively. These signals primarily represent the macroscopic activity of cerebral cortex. A single EEG signal is recorded by a single electrode, and is a spatially and temporally smoothed version of neuronally generated electrical potentials [23]. Intracranial EEG is invasive and difficult to place enough electrophysiological recording electrodes in a given cortical location to understand how the underlying neurophysiological processes give rise to the EEG. As a result, modelling is a way to bridge the EEG data and the underlying neurophysiological processes.

There are predominantly two approaches to build a mathematical model of the brain: the 'bottom-up approach' and the 'top-down approach' [24]. The bottom-up approach forms a macroscopic model of the brain or a brain region by coupling a network of microscopic models of individual neurons (e.g. Hodgkin–Huxley neurons). To obtain the spatial scale relevant to modelling EEG signals, the large number of neurons involved in this kind of model, as well as the heterogeneity and non-local connectivity of neurons, makes this largely impractical for understanding which key population-level neurophysiological variables influence macroscopic dynamics.

On the other hand, the top-down approach defines dynamical models based on phenomenological observations and emergent behaviour. This means that the mathematical model is built to model large-scale signals, such as EEG, through the interactions of neural populations, where the neural populations only coarsely model the average properties of the underlying detailed neuronal networks. Although it would be useful to model and infer the detailed microscopic neurophysiological processes, this is difficult to do as a result of the aforementioned experimental and observability issues as well as computational tractability. As a result, we focus here on inferring population-level neurophysiological processes using top-down models.

We focus on one of the most common NMMs used to model cerebral cortex, the Jansen–Rit model [9, 25]. The Jansen–Rit model was developed based on the model proposed by Lopes Da Silva et al [26], who proposed a two-population neural field model to study spontaneous EEG generation and alpha frequency rhythms, and reduced the complexity of the model by aggregating the activities within a single cortical column to perform analyses. Jansen et al [25] constructed a similar model of pyramidal neurons to show that the mechanisms that guide the evolution of evoked potentials and spontaneous EEG signal generation are the same. Jansen and Rit [9] extended the model with a third population of local excitatory interneurons that provide feedback to the pyramidal neuron population. The latter model has three populations: pyramidal neurons, inhibitory interneurons and excitatory interneurons.

The output of the Jansen–Rit model used in this paper is considered an approximation of the EEG/iEEG signal. The model is based on a mean-field approximation, which simplifies the complex behaviour of individual neurons by averaging their activity. Due to this, Jansen–Rit model does not provide the precise representation of the EEG signal. However, it is helpful in understanding the general behaviour of cortical neuronal populations, which provides useful insights of making prediction about the behaviour of the EEG signal under certain conditions. Thus, we apply the EEG signal to the output of the Jansen–Rit model.

Figure 1 illustrates the model. The spatially averaged post-synaptic potential of population $m$ arising as a result of input from pre-synaptic population $n$ is expressed as

$$\begin{equation}{v_{mn}}\left( t \right) = {\alpha _{mn}}\mathop {{{\mathop \int \nolimits }}}\limits_{ - \infty }^t {h_{mn}}\left( {t - t{^{^{\prime}}}} \right)\phi \left( {{v_n}\left( {t{^{^{\prime}}}} \right)} \right){\text{d}}t{^{^{\prime}}},\end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

where ${\alpha _{mn}}$ is the population averaged synaptic connectivity strength. The post-synaptic response kernel is given by

$$\begin{equation}{h_{mn}}\left( t \right) = \eta \left( t \right)\frac{t}{{{\tau _{mn}}}}{\text{exp}}\left( { - \frac{t}{{{\tau _{mn}}}}} \right),\end{equation} \tag{ 2 }$$

where $\eta \left( t \right)$ is the Heaviside step function and ${\tau _{mn}}$ is the population averaged synaptic time constant. Moreover, $\phi \left( v \right)$ is a sigmoid function that transforms the membrane potential to a firing rate given by

$$\begin{equation}\phi \left( v \right) = \frac{1}{2}\left( {{\text{erf}}\left( {\frac{{v - {v_0}}}{\varsigma }} \right) + 1} \right),\end{equation} \tag{ 3 }$$

where ${v_0}$ is a threshold and $\varsigma$ is the slope of the sigmoid function, which is also the variance of firing thresholds of the population assuming the firing thresholds follow a Gaussian distribution. The mean membrane potential of population $m$ is the sum of the synaptic contributions,

$$\begin{equation}{v_m}\left( t \right) = \mathop \sum \limits_{n = 1}^N {v_{mn}}\left( t \right).\end{equation} \tag{ 4 }$$

The convolution in equation (1) can be also written as two coupled, first-order, ordinary differential equations [27],

$$\begin{equation}\frac{{{\text{d}}{v_{mn}}}}{{{\text{d}}t}} = {z_{mn}}\phantom{\frac{{{\alpha _{mn}}}}{{{\tau _{mn}}}}{\phi _{mn}} - \frac{2}{{{\tau _{mn}}}} - \frac{1}{{\tau _{mn}^2}}{v_{mn}}}\end{equation} \tag{ 5 }$$

$$\begin{equation}\frac{{{\text{d}}{z_{mn}}}}{{{\text{d}}t}} = \frac{{{\alpha _{mn}}}}{{{\tau _{mn}}}}{\phi _{mn}} - \frac{2}{{{\tau _{mn}}}}{z_{mn}} - \frac{1}{{\tau _{mn}^2}}{v_{mn}},\end{equation} \tag{ 6 }$$

where ${\tau _{mn}}$ is a lumped time constant, and $m$ and $n$ indicate the pre-synaptic and post-synaptic neural population, respectively.

Furthermore, we can express the NMM in matrix-vector notation so the operations can be expressed more compactly. The NMM can be expressed in matrix notation as

$$\begin{equation}\dot x\left( t \right) = Ax\left( t \right) + B\vec \phi \left( {Cx\left( t \right)} \right)\end{equation} \tag{ 7 }$$

$$\begin{equation}y\left( t \right) = Hx\left( t \right) + v\left( t \right),\end{equation} \tag{ 8 }$$

where $H$ is the observation matrix, $v\left( t \right)\sim N\left( {0,R} \right)$ is the observation noise [10] and $y\left( t \right)$ is the membrane potential of the pyramidal population, which is considered to be the contributor to the generation of the EEG signal in our model. The matrix $A$ is defined by the membrane time constants,

$$\begin{equation}A\, = \,{\text{diag}}\left( {{{{\Psi }}_1} \ldots {{{\Psi }}_n}} \right),\end{equation} \tag{ 9 }$$

where

$$\begin{equation}{{{\Psi }}_n} = \left[ {\begin{array}{*{20}{c}} 0&1 \\ { - \frac{1}{{{\tau ^2}_{mn}}}}&{ - \frac{2}{{{\tau _{mn}}}}} \end{array}} \right].\end{equation} \tag{ 10 }$$

The matrix $B$ is the synaptic gains from internal inputs, which has the form

$$\begin{equation}B = {\text{diag}}\left( {0\;{{\text{a}}_1} \ldots 0\;{{\text{a}}_n}} \right).\end{equation} \tag{ 11 }$$

The matrix $C$ is the adjacency matrix that defines the connectivity structure of the model,

$$\begin{equation}C = \left[ {\begin{array}{*{20}{c}} 0&0& \ldots &0&0 \\ {{c_{2,1}}}&0&\;&{{c_{2,{N_x} - 1}}}&0 \\ \vdots &\,& \ddots &\,& \vdots \\ 0&0&\,&0&0 \\ {{C_{{N_x},1}}}&0& \cdots &{{C_{{N_x},{N_x} - 1}}}&0 \end{array}} \right],\end{equation} \tag{ 12 }$$

which is a matrix of zeros and ones, where one indicates there is a connection between two populations and zero indicates there is no connection.

In order to apply a nonlinear Kalman filter to this model to infer both the state vector and the parameters, we define a vector of parameters as $\theta = {\left[ {u\,{\alpha _{pe}}\,{\alpha _{pi}}\,{\alpha _{ip}}\,{\alpha _{ep}}} \right]^T}$. This set of parameters corresponds to the input $u$ to the model and the population averaged connection strengths between the three populations. Moreover, we assume the time constants of the model to be constant to simplify the estimation problem. The above parameters were combined with the state vector $x$ to form the augmented state vector,

$$\begin{equation}\xi = {\left[ {{X^T}{\theta ^T}} \right]^T}.\end{equation} \tag{ 13 }$$

The augmented state-space model is then

$$\begin{equation}{\xi _t} = {A_\theta }{\xi _{t - 1}} + {B_\theta }\phi \left( {{C_\theta }{\xi _{t - 1}}} \right) + {W_{t - 1}},\end{equation} \tag{ 14 }$$

where ${W_t}$ is Gaussian noise.

Before defining the novel LSTM-based filter for state vector and parameter estimation, this section outlines the Kalman filter applied as a benchmark. In particular, a filter referred to as the Analytical Kalman filter (AKF) is applied. This filter is highly stable and accurate and was developed in prior work [2, 10, 27] that evolved from deriving the Kalman filter for general nonlinear NMMs using the specific sigmoidal non-linearity in equation [15]. For the sake of brevity, we refer the reader to these prior works for greater mathematical insight. Moreover, links to code provided at the end of this paper give the exact specification of the implementation of the AKF. Here, a brief description of the main computations of the AKF are provided.

The aim of the AKF is to find the most likely posterior distribution of the augmented state given the previous measurements under Gaussian assumptions. Such a posterior distribution is characterised by its a posteriori state vector estimate and its covariance,

$$\begin{equation}\hat \xi _t^ + = E\left[ {{\xi _t}{ }|{ }{y_1},{y_2}, \ldots ,{y_t}} \right]{ }\end{equation} \tag{ 15 }$$

$$\begin{equation}\hat P_t^ + = E\left[ {\left( {{\xi _t} - \hat \xi _t^ + } \right){{\left( {{\xi _t} - \hat \xi _t^ + } \right)}^ \top }} \right].\end{equation} \tag{ 16 }$$

The AKF proceeds in two steps: prediction and update. During prediction, the prior distribution (obtained from the previous estimate) is propagated through the model equations. This step provides the so called a priori estimate, which is also a Gaussian distribution with mean and covariance,

$$\begin{equation}\hat \xi _t^ - = E\left[ {{\xi _t}{ }|{ }{y_1},{y_2}, \ldots ,{y_{t - 1}}} \right]\end{equation} \tag{ 17 }$$

$$\begin{equation}\hat P_t^ - = E\left[ {\left( {{\xi _t} - \hat \xi _t^ - } \right){{\left( {{\xi _t} - \hat \xi _t^ - } \right)}^ \top }} \right].\end{equation} \tag{ 18 }$$

During update, the a posteriori state vector estimate is determined by correcting the a priori state vector estimate with recorded (EEG) data by

$$\begin{equation}\hat \xi _t^ + = \hat \xi _t^ - + {K_t}\left( {{y_t} - H\hat \xi _t^ - } \right),\end{equation} \tag{ 19 }$$

where ${K_t}$ is the Kalman gain [18],

$$\begin{equation}{K_t} = \frac{{\hat P_t^ - {H^ \top }}}{{H\hat P_t^ - {H^ \top } + R}}.\end{equation} \tag{ 20 }$$

The a posteriori state vector estimate covariance is then updated by using the Kalman gain,

$$\begin{equation}\hat P_t^ + = \left( {I - {K_t}H} \right)\hat P_t^ - .\end{equation} \tag{ 21 }$$

After each time step, the a posteriori estimate becomes the prior distribution for the next time step, and the process repeats for each new recorded data point. The AKF requires the a postersiori state vector estimate and its covariance to be initialised at time $t = 0$. The special features of the AKF are that the state vector estimate and its covariance, along with the Kalman gain, were derived using the fully non-linear Jansen–Rit model.

While the AKF framework can be applied to multichannel EEG data in a multivariate fashion, for simplicity this paper focuses on estimating the state vector and parameters of the Jansen–Rit model using a single channel of EEG data (analogous to the output observation/measurement of the Jansen–Rit model) as input to the filter.

The RNN was introduced mainly to deal with time series data such as speech recognition and natural language processing [28]. An ordinary feed forward neural network is only processing independent data points. However, when data points depend on previous data point, in the case of time series data, an RNN incorporates the dependency by providing feedback to the next timestep. In this way, the final output depends not only by the input but also on the outputs of previous neurons.

LSTM is a variation of the RNN that has improved its performance [20]. Since gradient descent is applied to train neural networks, the gradient might explode or vanish after applying the sigmoid function over and over again, limiting to a certain number of discrete time steps to avoid this problem [29, 30]. However, the LSTM provides stability to bridge many more discrete time steps by enforcing constant error flow within special memory cells.

LSTM neural networks incorporate memory cells and gate units to convey useful information about previous states of the neural network to the current state. An LSTM layer consists of multiple recurrently connected memory blocks and three multiplicative gates, known as input, output and forget gates, that learn to open and close access to the error flow within the memory cell [20, 31]. For example, an input gate could use inputs from other memory cells to decide whether it should store information in its own memory cell.

An extension of the LSTM model is the bidirectional LSTM model that provide more accurate and improved results [32–34]. One of the shortcomings of RNNs is that they can only learn from the context of the previous time steps and not anything after the current time step. Bidirectional LSTM can process data in both directions with separated hidden layers, which are then fed to the same output layer [35]. Also, since the LSTM is free to access as much or as little of the data within the given time window as needed, predefining a specific time window for the model is not required [31].

Here, a bidirectional LSTM filter is constructed that takes as input simulated or real EEG and predicts the state vector and parameters of the Jansen–Rit model as well as simulated or real EEG. Similar to the AKF, while the framework that follows could potentially be applied to multichannel EEG data in a multivariate fashion, this paper, for simplicity, focuses on using the LSTM filter to estimate the state vector and parameters of the Jansen–Rit model for a single channel of EEG data.

The structure of the LSTM model is designed to achieve high performance while maintaining good time efficiency. Since the only non-linearity of the NMM is the firing rate function in equation [3], and to achieve high performance while predicting parameters given the observation, we would expect that at least two layers of artificial neural network are required. A grid search experiment was conducted to test the number of neurons in each layer. The selected structure was 128 and 32 neurons in layers 1 and 2, respectively.

As real EEG data provides no ground truth about the neurophysiological variables that correspond to the state vector and parameters of the Jansen–Rit model, the LSTM filter is trained on simulation data was required. This allows the prediction of the state vector and parameters of the model when inputting either simulated or real EEG into the LSTM filter. It is important to generate a variety of training data with different patterns and a wide range of parameters so the model can learn the relationship between signal patterns and how the state and parameters change. Furthermore, it is also crucial to generate data containing oscillations, since there are likely to be many solutions for a constant signal, or fixed point. This is because when there is no oscillation, different combinations of changes in either the external input, input from the inhibitory or excitatory interneurons may result in the same observation signal. However, changes of the parameters could affect the ranges of external input that generate oscillations, which makes it difficult to find the desired range of external input when there can be different combinations of parameters. This is because the time constants affect the transformation from firing rate to post-synaptic membrane potential thus affecting the ranges over which the population produces oscillations. Thus, an effective method to generate oscillatory data with a wide range of parameters is required.

Different combinations of inhibitory and excitatory time constants can generate data with different waveforms, frequencies and amplitudes. Normally, the ranges of both excitatory and inhibitory time constants are considered between 10–60 ms [36]. In addition, external input also plays an important role in data generation. External input to the model is also affected by the time constant, such that the range of the external input that can produce oscillations is affected by different combinations of time constants, but the exact range for each combination is unknown. To determine the range of external input and whether the given combination of time constants is able to generate oscillation, an automatic data generation process is designed.

Statistical hypothesis testing is used to determine whether an oscillation exists in a given time-series recording. Two tests are used. If the simulated EEG signal is noise-like, then it will have a Gaussian distribution so the Anderson–Darling test [37] is used to test if the data is not Gaussian distributed, which means there could be an oscillation. The other test is the Ljung–Box test [38], which tests whether any autocorrelations of the time series are different from zero. Either of the tests can help determine whether the generated data contains an oscillation and, since it is more important to ensure there are no false positives (data is Gaussian distributed or data has non-zero correlation, but we reject this hypothesis) to avoid data being generated without including oscillations. We set the significance level to $\alpha = {10^{ - 4}}$.

A flow chart used to generate the training data is shown in figure 2. We generate a short recording first to check whether this combination of time constants and external input is able to generate data containing oscillations. If not, then the external input is increased until it does. If it still does not produce oscillations for 15 consecutive increases of the input, the current combination of time constants is considered to not be able to produce the desired data and so the next combination is used. If the generated data includes oscillations then the upper and lower bounds of the external input are recorded, and training data is generated with the recorded parameters.

Zoom In Zoom Out Reset image size

The size of the dataset is as large as possible within the size of physical memory, which means that the gap between each combination is kept relatively small and as many examples as possible are provided for training [39]. The generated data was shuffled randomly and divided into training, validation and testing data with ratio of 80:10:10.

To improve the performance of the model with different amplitudes, the dataset is standardised. As in many cases with either real or simulated data, the raw amplitude of the EEG data may be on arbitrary scales. To avoid potential issues due to variations in amplitude scale or values out of the training range, The mean and standard error were computed for each feature and target variable, and standardisation is implemented for each of them using $\frac{{x - \mu }}{\sigma }$. The other reason for doing this is to ensure the loss is reduced for all variables without any bias, as different amplitudes might cause the loss of some of the target variables to be weighted more than others. For example, if one variable has a higher amplitude than the other, its error is higher and thus will be reduced more.

The training of an artificial neural network aims to minimise the loss function over all training data. By default, the loss function is the mean squared error of all target variables, which are the state variables and parameters in our case. However, it is not enough for the model error to be minimised only on the training data, as we would like the LSTM model to learn whether the predicted state variables and parameters can represent the observations. More importantly, it should also learn to behave like the NMM. Given an observation signal, the trained LSTM model should be able to predict the state vector and parameters that can recreate the observation signal using the NMM. Since the generation of the training data is based on setting the discrete parameters, even though we can lower the step within each parameter, it is still possible that the model cannot learn the mathematical relationship between all state variables and parameters. Thus, we link the LSTM model with the NMM.

To link the LSTM model with the NMM, we customise the loss function to allow the LSTM model to know the mathematical relationship between the observation and the state. Since the LSTM model is able to produce the state for the next timestep $t$ given the state of the current time step $t - 1$, it is possible to generate all $t$ states given the state at the current timestep. By comparing the state at time $t$ with the state predicted by the LSTM model, we can know the error between the LSTM model prediction and the NMM prediction. If we can link the NMM to the loss function, the LSTM model can learn to minimise both the error of the LSTM predicted values and the NMM predicted values

The squared error at a single timestep is defined as

$$\begin{equation}{\text{Squared Error}} = \left[ {{{\left( {{\xi _{\text{t}}} - {{\hat \xi }_t}} \right)}^2},\,{{\left( {{y_t} - H{{\hat \xi }_t}} \right)}^2}} \right]\end{equation} \tag{ 22 }$$

where $\xi$ is the augmented state, $y$ is the observation of the training data and $\hat \xi$ and $H{\hat \xi _t}$ are the augmented state and observation predicted by the LSTM model, respectively. The Squared Error term compares the predictions with the true values. By minimising the error, the predictions of both states and measurements will be closer to the truth.

We then add to the loss function the model error, which is the error between the state of the LSTM model and the state generated by the NMM [40, 41],

$$\begin{equation}{\text{Model Error}} = \left[ {{{\left( {{{\hat \xi }_t}\! -\! \left( {{A_\theta }{{\hat \xi }_{t\! -\! 1}} \!+\! {B_\theta }\phi \left( {{C_\theta }{{\hat \xi }_{t \!-\! 1}}} \right)} \right)} \right)}^2},{\text{ }}0} \right].\end{equation} \tag{ 23 }$$

With the augmented state-space model, the estimated augmented state vector at the current time step ${\hat \xi _{t - 1}}$ can be converted to the next time step. Minimising the squared error between ${\hat \xi _t}$ and the converted state vector ${A_\theta }{\hat \xi _{t - 1}} + {B_\theta }\phi \left( {{C_\theta }{{\hat \xi }_{t - 1}}} \right)$ links the estimation to the NMM, since the training has to minimise the error between its own prediction and the state generated by the NMM.

Note that, although ${B_\theta }$, ${C_\theta }$ and $\phi$ are considered to be constant, the matrix $A$ can vary depending on the time constants; i.e. the $A$ has to change over time. In the customised loss function, $A$ is calculated at each timestep so the change in state and parameters can be reflected over time.

The last thing aspect of the loss function is the rate of change of the connectivity strength parameters, as they are considered to be slowly changing parameters compared to the membrane potentials and simulated EEG. It is possible to add the standard deviation of the parameters to the loss function to limit the amounts that they change. However, the parameters need to be adjusted rapidly when they differ substantially from the NMM. Thus, the standard deviation is combined with the model error as

$$\begin{equation}{\text{Std = }}\;s\left( \alpha \right)\left[ {{{\left( {{{\hat \xi }_t} - \left( {{A_\theta }{{\hat \xi }_{t - 1}} + {B_\theta }\phi \left( {{C_\theta }{{\hat \xi }_{t - 1}}} \right)} \right)} \right)}^2},\,0} \right]*k\end{equation} \tag{ 24 }$$

$$\begin{equation}s\left( \alpha \right) = {\left[ {0, \ldots 0,{\text{std}}\left( \alpha \right),0} \right]^T},\end{equation} \tag{ 25 }$$

where $\alpha$ is the vector of the four connectivity strengths parameters and $k$ is an adjustable weight that is set to 0.1 as default.

The final loss is the summation of equations (21)–(23).

On top of the testing data split from the training dataset, another testing dataset with the same range of parameters was generated, but the steps of ${\tau _e}$ and ${\tau _i}$ are set to be smaller. This dataset avoids using the exact parameters of the original training dataset.

We tested both the AKF and the LSTM model on the same dataset. To evaluate the impact of incorrect initialisation of the AKF, two settings were used. The first setting was perfectly initialised Kalman filters that were initialised to be exactly the same as the parameter initial values. The second setting was default parameters, where the excitatory and inhibitory time constants were 0.02 s and 0.01 s, respectively, which corresponds to the parameter values that generate a strong alpha-like (8-12 Hz) rhythm [9].

To ensure all state vectors and parameters are on the same amplitude scale numerically, the testing results are standardised. We compute the root mean squared error (RMSE) between the truth and the prediction for each state variable and parameter, and find the difference between the AKF and the LSTM model result. The RMSE for each parameter is

$$\begin{equation}{\text{BCM}}\;{\text{ = }}\;\sqrt {\mathop \sum \limits_{i = 1}^n \left.{{\left( {\frac{{{x_i} - \mu }}{\sigma } - \frac{{{{\hat x}_i} - \mu }}{\sigma }} \right)}^2}\right/n} ,\end{equation} \tag{ 26 }$$

where $\mu$ and $\sigma$ are the mean and standard deviation of the testing dataset for each parameter, respectively, and ${x_i}$ and ${\hat x_i}$ are the truth and prediction, respectively. Equation [26] is used for each estimation method to extract the raw RMSE and the values are compared across methods.

Testing the robustness of the model is also of great importance, since typically there will be strong noise in real EEG data. As the standard deviation of the testing observation data is around 40, we add Gaussian noise with a mean of 0 and a standard deviation of 4, representing 10% Gaussian noise added to the testing data. Both the LSTM model and the perfectly initialised Kalman Filter are tested using the noisy testing data.

We also test the LSTM model on randomly generated data involving time-dependent parameter changes. The time constants are randomly chosen from the range 0.01–0.06 s, and the model is simulated with the selected time contents held fixed for 5 s before the next combination of time constants. Note that there is a 5 s buffer between each 5 s fixed time constant simulation segment, where we connect the parameter values between segments with a straight line to simulate transitions between different brain states, where the parameters are not constant. We assume the parameters are slowly-changing, so a straight line between different stable states would be enough for simulation purposes. Since the training was done on constant parameters with noise only, testing on randomly generated data is needed to know whether the LSTM filter can work in a scenario where parameters are time varying.

Apart from simulation data, we also use intracranial EEG (iEEG) data collected from a patient with epilepsy chosen from data recorded continuously from 15 patients for up to three years [42]. The patients were implanted with 16 intracranial electrodes around their seizure onset areas with sampling rate 400 Hz. The electrode array was wirelessly connected to a portable device. Seizures were automatically detected by the advisory device and were reviewed by clinical experts. We denote the iEEG data from this study as 'real data'.

To demonstrate the practical utility of the LSTM filter, it is applied to a 1 h intracranial EEG recording containing an epileptic seizure to show whether or not we can observe a transition between normal and ictal brain states. Although the actual values of the state variables and parameters are not known, there is an obvious brain state change with seizure onset, so we can see if the LSTM model may be used to detect or even predict a seizure. Since normally iEEG data becomes more rhythmic when it enters a seizure compared to normal activity, we expect the parameter estimates to become stable after entering the seizure state.