Computational Sensitivity Analysis on a Mathematical Model of Epileptic Seizures

Temporal lobe epilepsy is a serious neurological disorder characterized by complex partial seizures, which are thought to originate in the hippocampus. Ordinary differential equation modeling has been used to describe changes in membrane potential of excitatory and inhibitory cells in order to gain insight into seizure propagation. In the current study, a system of ordinary differential equations based on previous modeling is used with distinct biologically reasonable values for membrane capacitance in order to determine model sensitivity to that parameter. Because delay differential equations are used in the model, sensitivity is investigated computationally by examining the variation in output relative to various inputs. Membrane capacitance was found to affect model predictions and whether groups of cells exhibited the same behavior after a certain period of time. Hence, membrane capacitance is a critical parameter when modeling changes in membrane potential and should be incorporated clearly. Changes in model output as a result of changes in a time delay parameter are also investigated. Keywords-epilepsy, membrane potential, sensitivity analysis, mathematical modeling


I. INTRODUCTION
Temporal lobe epilepsy is a serious neurological disorder characterized by complex partial seizures, which are thought to originate in the hippocampus [1], [2], [3], [4], [5].The mechanism of the propagation of these seizures is not completely understood, and there is no single pathology even within temporal lobe epilepsy [2].Surgical treatments are often recommended and beneficial for patients with epilepsy [5], but the success of these treatments is dependent on the completeness of the tissue removal and the surgical approach [6], [7], [8], [9] as well as a strong understanding of a patient's clinical history [10].At this time, surgical techniques may be the most effective treatment strategies from both cost and outcome perspectives [11], [12].Surgical techniques carry inherent dangers and mathematical modeling may aid in developing alternative treatment strategies with positive outcomes [13], [14].
Often in modeling, fixed parameters are based on actual measurements.One parameter used in modeling based on the Hodgkin and Huxley model [34] is membrane capacitance.Most plasma membranes have a capacitance around 1 µF/cm 2 [35] which suggests that a value of C = 1 is reasonable.Similar models have explicitly used C = 1 [23], but there are others that use C = 1 but do not explicitly include the parameter [24], [25].However, a value close to (but not equal to) 1 may result in different model output than a value of exactly 1.The study in Gentet et al. [36] looked specifically at membrane capacitance in different neurons.All of their values were indeed around 1, but those values for hippocampal neurons (from eight Wistar rats) were 0.92±0.08µF/cm 2 .The other neuron mean values ranged from 0.85 to 1.06 µF/cm 2 .Additionally, the measured value of membrane capacitance in neurons may depend greatly on the measurement techniques used [37].Hence, a local sensitivity analysis on the model in the current study is performed to ascertain the effect of changes in C on the model output.If the model output changes greatly with small changes in the value of C, this parameter is critical to include when making model predictions or modifying the equations.
In the current study, a previously developed model from Larter et al. [25] is modified to be more defined with respect to units and to be more consistent with previous literature.Additionally, the revised model incorporates a more specific time delay.In the course of modifying the model from Larter et al. [25], a parameter denoting membrane capacitance is reintroduced in the model.This parameter, referred to as C, is expected to have a value close to 1 µF/cm 2 .A sensitivity analysis is performed in order to determine if values of C that are close to but not equal to one will have a large impact on model predictions.Because the model involves delay differential equations, the sensitivity investigations will be performed computationally.The impact of small changes in the value of C are important since the Larter model essentially eliminates the parameter C by setting it equal to one [25].Hence, this study is intended to model membrane potential changes in the hippocampus using foundational equations [34], [38] in context with the addition of inhibitory neurons in a similar way to the study in Larter et al. [25].Further, this study investigates the importance on the accuracy of the value of the membrane capacitance parameter when making predictions using this model.The parameter for time delay is also varied in the process of the investigation.

A. Model Structure and Development
The study in Larter et al. [25] used a system of three ordinary differential equations to describe membrane potentials of excitatory and inhibitory neurons in each node of a 36×36 lattice.The lattice was intended to represent the CA3 region of the hippocampus because this region is a common location of the epileptic focus.The model in Larter et al. [25] is based on equations describing voltage oscillations in the barnacle giant muscle fiber [38] which itself is based on a model by Hodgkin and Huxley [34].The foundation for the pyramidal (excitatory) nerve equations in Larter et al. [25] is a model of a single voltage-dependent conductance in which the calcium (Ca 2+ ) system is assumed to be much faster than the potassium (K + ) system as described in the study by Morris and Lecar [38].Larter et al. [25] adds one additional differential equation in each node for the inclusion of inhibitory neurons.
The units of the model in Larter et al. [25] are primarily arbitrary, but the units in the source equations are specifically defined.The models in earlier studies [34], [38] used a particular parameter, C (membrane capacitance, µF/cm 2 ), but the model in Larter et al. [25] essentially eliminates this parameter by setting C = 1, which is a reasonable value [35], [36] and has been used in other similar models [23].Additionally, Larter et al. [25] did not use delay differential equations but did incorporate a time difference when modeling K + transfer from surrounding nodes.
The model in the current study is similar in concept to the model in Larter et al. [25] in that it incorporates the effect of inhibitory neurons, but the current model uses two equations per node for this effect instead of one.These equations are changed in order to model the membrane potential of inhibitory cells the same way the potential of the excitatory neurons are described in Larter et al. [25] and to be more similar to previous models [23], [34], [38].Sodium channels, Na + , are represented (instead of calcium channels) as in [34].Additionally, a continuous delay is used when describing K + transfer from surrounding nodes, and membrane capacitance (C) is reintroduced as an explicit parameter.Both of these changes were made to make the model more physically representative.The updated equations for each The states at each node i are V i , the membrane potential of the pyramidal cells; Z i , the membrane potential of the inhibitory cells; and W i and N i , relaxation factors.Descriptions of the parameters (other than C) and their fixed values are presented in Table I.
Biomath 4 (2015), 1508141, http://dx.doi.org/10.11145/j.biomath.2015.08.141 Equations ( 1) and ( 2) are the same as in the Larter model [25] with the parameter for membrane capacitance, C (µF/cm 2 ) explicitly included and a temperature scaling factor removed.Equation ( 1) is directly from the model from Morris and Lecar [38] with an additional loss term, α inh (Z i )Z i , representing the effect of the inhibitory neurons on the membrane potential of the pyramidal cells.The equation for Z i from Larter et al. [25] was replaced in this study with a differential equation (equation ( 3) above) that is similar in structure to equation (1).In equation ( 3), the excitatory effect on the inhibitory cells is modeled through an additive term, α exc (V i )V i , which is similar to the way that effect was modeled by Larter et al. [25].The rest of equation ( 3) is directly from the Morris and Lecar model [38].Both equations ( 1) and ( 3) are based on the Morris and Lecar model of single voltage-dependent conductance in which the Ca 2+ system is assumed to be much faster than the K + system [38], and a similar assumption is made here for the Na + channels.Equation ( 4) is similar to equation ( 2) and is added in order to model a relaxation factor for the inhibitory neurons in the same way a similar relaxation factor is modeled for the excitatory neurons.
The additional functions used in equations ( 1)-( 4) are defined below: where X i is either V i or Z i .The functions defined in equations ( 5)- (10) were used the Larter model [25] as taken from the model by Morris and Lecar [38].Other studies have used different functions that are also sigmoidal [23], [24], [27].The notation of the constants in equations ( 5)-( 10) has been changed from the Larter model [25] in order to avoid confusion.Lowercase letters are now used so that constants will not be confused with state values at particular nodes.The parameter values used in equations ( 5)-( 10) are presented in Table II.The additional functions added in the current study, τ n (Z i ) and n ∞ (Z i ), are defined identically to functions with a similar purpose.Equations ( 5) and ( 6) represent the fraction of open voltage-regulated Na + channels, and the effect of excitation or inhibition is modeled in a similar way in equations ( 8) and (9).Equations ( 5), (6), and (10) were also used in the Larter model [25] and taken directly from the study by Morris and Lecar [38].More information on the functions used in equations ( 5)-( 10) can be found in the study by Larter et al. [25].
The model in Larter et al. [25] did not use delay differential equations but did incorporate a time step when modeling the change in membrane potential due to extracellular potassium.The quantity V K i in the current study and in Larter et al. [25] models how the membrane potential of surrounding nodes affects node i.In the current study, however, this effect is modeled using a delay term as described below where j indicates the number of a node adjacent to node i and X represents either V or Z.Each node is assumed to contribute equally to an adjacent node and periodic boundary conditions are used.Hence, X j (t − τ )/6 is the delayed contribution of node j to node i for either excitatory (V ) or inhibitory (Z) neurons.The lattice used in the current study is similar in structure to that used in the Larter et al. study [25] in that it is composed of hexagonal "nodes."The current study uses a smaller lattice as shown in Figure 1.Nodes are numbered beginning with the upper left corner and proceeding left to right and on to additional rows.This illustration also includes a closer view of the six nodes surrounding any internal node.Conceptually, the lattice in the current study does differ from that of other studies [23] in that the inhibitory cells are considered to affect other inhibitory nodes surrounding them.The lattice is intended to represent populations of cells (as with [23]) which conforms to the strategy of using variables not at a microscopic but at a macroscopic level as used in the study by Kaneko [39].
Although the authors from Larter et al. [25] did not use τ , their model used an integration time period, t int , which functioned in a similar capacity.Larter et al. discovered that the lattice became synchronized, i.e., the difference between node values was negligible, if t int fell below a particular threshold (13.5 in arbitrary units) and used various values of t int for their simulations (including a value as high as t int = 48).Because the current study uses a somewhat different model framework for this time value and because τ is specifically defined in seconds in the current model, τ = 1, 2, 5, and 10 s were used for the The leak conductance parameter V L was altered in order to make the solutions exhibit more cyclic behavior, i.e., so that the node potentials did not continue to increase or to decrease.

B. Computational Methods
Due to computational restrictions, a smaller lattice size (20 × 20) was used than in the model from Larter et al. [25].The differential equations representing the system, (1)-( 4) were solved numerically for C = 0.8, C = 1.0, and C = 1.2 and various values of τ .All simulations were performed using Wolfram Mathematica, version 9.0.The state values for t < 0 were set equal to their respective initial conditions, and initial conditions did not vary across the lattice.The simulations were run for 350 seconds in order to allow the curves to approach a repeated cycle.The last 150 seconds were used for plotting and are intended to represent a general time period of that length.
The Larter model on which the current study is based investigated the difference in membrane potential between nodes as a measure of synchronization of the lattice [25] because synchronization may be an indication of epileptic phenomenon [1], [2], [32], [27].A similar analysis was conducted in the current study, and nodes 43 and 358 were chosen for investigation because they were not on the boundary and relatively distant from each other as is shown in Figure 1.

III. RESULTS
In order to determine how the reduction in lattice size affected results, computational solutions were generated using a 6 × 6, 10 × 10, and 20 × 20 lattice and compared.Values of C = 1, and τ = 1 were fixed for these solutions.Nodes 9 and 28, 23 and 78, and 43 and 358 were investigated for the 6 × 6, 10 × 10, and 20 × 20 lattices, respectively.Figure 2 contains graphs of the differences in the membrane potentials between the investigated nodes, and differences can be seen in the node differences for both excitatory and inhibitory potentials.
Selected plots of the computational solutions to equations ( 1)-( 4) are presented in Figures 3-9 An examination of the results for the membrane potential in excitatory neurons (presented in Figure 3) would suggest that changes in the membrane capacitance (C) result in changes in the amplitude and phase of the oscillations.Also, it appears that changes in the delay parameter (τ ) lead to a change in the period/frequency of the oscillations in the membrane potential, as well as the amplitude.Variation in the membrance capacitance and the delay parameter lead to similar changes in the membrane potential in the inhibitory neurons (presented in Figure 4), although changes in the capacitance also seem to affect the baseline of the oscillations.
Differences between predictions at node 43 and node 358 are presented in Figure 5.These differences are plotted in order to ascertain if the two nodes are predicted to exhibit the same behavior with time.The magnitude of these differences is similar in scale to model output for excitatory and inhibitory cells, which indicates that the model is predicting different values at the investigated nodes when time is fixed.
The membrane potentials for the inhibitory cells shown in Figure 4 and for the pyramidal cells shown in Figure 3  these figures represents the behavior of inhibitory cells versus the behavior of excitatory cells at an individual node and specific value of membrane capacitance at a particular value of τ .The shapes of these curves remained fairly consistent with changes in the value of C but the cycles were more erratic for increasing values of τ .Differences in predictions for inhibitory nodes and excitatory nodes are plotted in Figure 10, and the highest differences (in magnitude) between nodes are presented in Table III for each investigated value of C and τ .

IV. CONCLUSIONS
A model describing the change in membrane potential in excitatory and inhibitory neurons in the CA3 region of the hippocampus [25] was modified to more explicitly contain a parameter representing membrane capacitance.The revised model uses ordinary differential equations for both excitatory cells (V i ) and for inhibitory cells (Z i ), and revisions included adapting equations for Z i to be more similar to those for V i .The revised model also incorporates delay differential equations to allow for signaling between nodes.A smaller lattice size was used in the current investigation than in the previous model [25], and variation between nodes does appear to increase with lattice size as shown in Figure 2.
In order to ascertain the importance of an explicitly defined capacitance parameter, model output was generated for various values of C and τ .As shown in Figures 3-10, values of C above and below 1 result in different model output for values of V i and Z i and affect the level of synchronization of the lattice the model predicts.
The model output is sensitive to changes in the value of C (for various values of τ ) which indicates that the accuracy of this parameter does affect model predictions.Hence, these results show how critical the choice of C is because C = 0.8 may be more realistic a value for C for hippocampal neurons than the more commonly used C = 1 [36].Model output is also affected by values of τ , with larger values of τ resulting in more oscillation.
Assumptions were made in the process of generating simulations, such as choosing similar parameters for inhibitory and excitatory nodes.Parameters such as channel conductances may vary and may affect model output [24].In this study, all parameters were fixed except for C and τ in the local, visual sensitivity analysis.Because the sensitivity analysis was local and visual, changes in other fixed parameters may affect output results.However, changes in model predictions were seen for different values of C which shows that for some sets of parameters, the model does depend on the value of C, and hence, the results indicate that C is important to accurate model output.The model also appears to be sensitive to the parameter τ and additional investigation on that parameter may be necessary as well.
Another aspect of the current study that may affect simulation output is the size of the lattice.The lattice used in the current study had fewer nodes than in the Larter model [25].Because the nodes represent groups of cells rather than individual neurons, the 20 × 20 lattice can be viewed as a coarser representation of the same system.A reduced lattice was also used in the current study for simplicity because the introduction of the delay and modeling inhibitory node communication increased computational difficulty.
The model in the current study describes changes in membrane potential in excitatory and inhibitory cells in the hippocampus and is based on previous models [25], [34], [38].This system of differential equations is intended to help understand epileptic seizures, and innovative treatment strategies may be developed through this understanding.The model is rooted in the specifics of the Hodgkin and Huxley model [34] and later simplifications of that model [38], but information may be learned about the overall system by using these simplifications [14].Additionally, the number of elderly patients with epilepsy and seizure disorders is likely increasing but access to relevant tissue samples can be problematic [41].Predictive mathematical models based on physiological mechanisms are potentially useful in identifying changes with age that affect seizures and hence may be useful in identifying treatment strategies for epileptic patients of various ages.The structure of the model may not reveal all aspects of epilepsy at this point, but more can be learned through refining and further investigation of this system through mathematics.

Biomath 4 (Fig. 1 .
Fig. 1.Illustration of the lattice used in the simulations. . The solution curves correspond to membrane potentials and their sensitivities for nodes 43 and 358 on a 20 × 20 lattice.The membrane potentials of excitatory cells at nodes 43 and 358 versus time are presented in Figure 3. (The choice of nodes 43 and 358 for investigation is discussed in Section II-B.)Similarly, the membrane potentials of inhibitory nodes versus time are presented in Figure 4. Three curves are presented in each of these figures which represent solutions using three different values of C. Different solution curves are apparent for the three values of C, and the largest range of excitatory output is present when C = 0.8 for most of the investigated values of τ .Additionally the largest magnitude of inhibitory output in the investigation occurs when C = 0.8 for most investigated values of τ .For cases (such as τ = 5) when the output ranges are larger for simulations using other values of C, there are visible differences in model output for various values of C.

TABLE II DESCRIPTIONS
AND VALUES OF PARAMETERS USED IN EQUATIONS (5)-(11).ALL VALUES BELOW WERE USED IN THE MODEL IN LARTER et al [25] ALTHOUGH THE NOTATION HAS CHANGED IN SOME CASES.VARIOUS VALUES WERE USED FOR SOME PARAMETERS IN LARTER et al. [25], BUT ONLY ONE WAS SELECTED FOR THIS INVESTIGATION.THESE VALUES ARE NOTED WITH AN ASTERIX (*).*simulations in order to visually evaluate sensitivity.