User:Jrm228/Final Term Paper

From BIOL 300 Wiki
Jump to: navigation, search

Author: James McGinnity

Introduction: Modeling Excitatory and Inhibitory Neurons

The Problem

In the past few decades, the combination of massive leaps in technology and of an increase in the accessibility of scientific data has led to incredible strides in modeling neuronal networks. However, the complexity of pre-existing models such as the Hodgkin-Huxley model has made them computationally inefficient and prevented a more complete understanding of the brain from being developed [1][2][3]. One researcher, Dr. Eugene Izhikevich, created a computationally efficient system that would still be biologically plausible by reducing the number of parameters and state variables to only the most influential parts of the system [4]. This model can be run in real time for tens of thousands of spiking neurons, but there are a variety of differences between this model and conventional models, and so there is a concern that Izhikevich’s model is not accurate enough to provide realistic simulations of brain activity.

Background and Hypothesis

There are two parts to demonstrating the accuracy of Izhikevich’s model that this paper will explore. The first part will involve reproducing the figures that he generated in his paper, which show the changes in the model for several different types of firing in neurons [4]. This will be compared to a Mathematica model of the Hodgkin-Huxley set of equations that is demonstrated in chapter 7 of our textbook [5]. The second part will expand upon both models by way of a bifurcation analysis, which can demonstrate the degree of difference between the two models to a much greater extent [6] The model proposed by Izhikevich does not lose a significant amount of data when compared to the Hodgkin-Huxley model, and the bifurcation analysis will demonstrate that the results of both models are similar.

Summary of Results

Our models, shown in the next section, were fairly similar on the surface. Both showed the spike at the same time, and had the same post-spike dip, but the Hodgkin-Huxley model had a lot more little bumps and irregularities that are characteristic of a more accurate model. We were unable to make a significant conclusion on our hypothesis, so we moved to conduct our extension through a bifurcation analysis. We modeled the nullclines and equilibrium points of the membrane potential of both sets of equation. This analysis ultimately demonstrated that a substantial difference existed between the two models.

Model Description

State Variables and Parameters [4]

State Variables

The following table lists the two state variables in the model. The initial values are for the first neuron, and are subject to change for every other neuron.

Variable Description Initial Value
v membrane potential -65 mV
u membrane recovery variable 0


The table below defines the five parameters present in the model. The first two, a and b, will influence the initial values of the state variables, and, together with the last, I, will influence the first derivatives. The other two, c and d, will effect the reset value when the spike occurs.

Variable Description Default Value RS IB CH FS LTS RZ
a timescale of membrane recovery variable 0.02 0.02 0.02 0.02 0.1 0.02 0.1
b sensitivity of recovery variable to fluctuations in membrane potential 0.2 0.2 0.2 0.2 0.2 0.25 0.25
c after-spike reset value of potential -65 mV -65 mV -55 mV -50 mV -65 mV -65 mV -65 mV
d after-spike reset value of recovery 2 8 4 2 2 2 2
I injected or synaptic current 10 10 10 10 10 10 10
  • Note: RS: Regular Spiking; IB: Intrinsically bursting; CH: Chattering; FS: Fast Spiking; LTS: Low Threshold Spiking; RZ: Resonator.

Assumptions [4]

There are several key assumptions that must be made in order to model the neuron with this set of variables and equations, as well as to make the extension;

  1. With regards to all five parameters, it is assumed that they accurately represent the cortical neuron as it would appear in a natural state.
  2. The models can be reasonably compared by the membrane potential. Both models have other state variables, but they are not consistent between the models, which makes them impossible to compare.
  3. A bifurcation analysis of regular spiking behavior is sufficient in demonstrating the properties of both models. The Hodgkin-Huxley changes the state variables substantially between different types of behavior, but the paper model only alters parameters.

Inputs to Model

The following are two equations describing the system and two conditions that are imposed upon the equations. After each, there is a description of what each component means in the context of the mathematical model[4].

Equation 1:

 \frac{dv}{dt} = 0.04v^{2} + 5v + 140 - u + I

This equation is representative of the membrane potential of the cortical neuron, which is the primary tool of sending out the waves we are interested in. The first three components are grouped together as a quadratic equation for the membrane potential itself, indicating that as the membrane potential grows, it speeds its rate of growth. The fourth term, u, describes the recovery variable that is shown in equation 2. It limits the rate of change of the membrane potential slightly, but not enough to stop it. Finally, the term I is representative of the injected or synaptic current in the cortical neuron. If there is a greater current present, then the potential will increase as well.

Equation 2:

 \frac{du}{dt} = a ( bv - u )

The membrane recovery variable is shown above as u, and is determined by this equation. The terms a and b are the parameters presented in the Parameters section above, and determine the rate at which the recovery variable slows growth. As u increases, the rate of change decreases, but the membrane potential increases so much more rapidly that the recovery variable continues to grow over time.

Condition 1:

If  v \ge 30 , then  v = c

This condition allows for a reset after each spike so that the signal can be sent out again. The parameter c represents the reset value for the membrane potential to return to after each spike, and the reset is triggered when the potential goes above 30 mV.

Condition 2:

If  v \ge 30 , then  u = u + d

At the same time, the value for the recovery variable must be changed. However, while the change to the membrane potential in the other condition always vastly decreases the value of v, this condition increases the value of u so that it can aid in swiftly returning v to its original, pre-spike level.

Simulation Method

In order to conduct our model, we had to first demonstrate that a simpler version of the Mathematica code would run correctly. We started by creating a model that would plot the first graph of Figure 2 of our paper by using the above equations and conditions with the parameter values of a Regularly Spiking neuron. This used an NDSolve to determine the equations for v and u, which we then plotted.

After confirming this, we then set out to create a similar model from Chapter 7 through use of the Hodgkin-Huxley equations in order to compare our model to a biologically accurate model. We used the code from the text, then altered the parameters until it approximated our model results.

Finally, in order to test our hypothesis and extend the original paper's model, we performed a bifurcation analysis on both sets of data. For the paper model, we used a Solve with the two equations in a quadratic form and without the WhenEvents that allowed for our conditional spike event. Chapter 7 was slightly trickier, as it was a four-dimensional system, but after applying a Solve and using our altered parameters, we were able to create nullclines and equilibrium points for it as well. We compared the two sets of results to each other to determine the validity of our hypothesis.

Basic Model of Spiking Neurons

We will go over the two different models for spiking neurons, with special attention given to the discrepancies between our model and the model that Izhikevich published in his paper. While the models have enough differences to comment on, the more extensive comparison will come with the bifurcation analysis following that section. For more information on how the figures were created, please consult the Mathematica notebook file at the following link; Basic Model

Izhikevich Model (Original Model)

The paper we are examining had several different graphs, but the one that drew the most attention was a graph with eight different types of neuron firing demonstrated by his code. This graph is shown below;

Figure 1: The graph comes from the paper by Izhikevich, and shows the known types of neurons corresponding to different values of the parameters. This is the figure that will be replicated later on, and, as with all of the figures, the discrepancies will be addressed outside of the captions.

By displaying the effects that each parameter has in sync with each other, Izhikevich is able to demonstrate the full scope of the model that he claims is not only much more computationally efficient than the Hodgkin-Huxley model, but also just as biologically plausible. The hypothesis that we are testing is whether or not that claim is accurate, and that means we need to replicate Figure 1.

The first step to replicating Figure 1 is to make a representation of the equations for a single set of parameters. To this end, we used NDSolve to input the parameters and apply them to his equations, then evaluated the results and graphed them to produce Figure 2.

Figure 2: This graph is of a single, regularly-spiking cortical neuron based on the code provided by Izhikevich in his paper. The value of the membrane recover variable, u[t], is shown above the membrane potential, v[t], though this will not be shown in the final replication.

There are several discrepancies between his graphs and ours that should be addressed. First, in this representation, the membrane recovery variable u[t] is displayed alongside the membrane potential v[t]. This was done to represent the reset of the system after each spike, which is clearly shown with u[t]. In addition, there is only one spike here compared to several spikes in the original paper model. This was a result of the time scope that is examined in each. His model used 0.1 milliseconds, while ours only went for 0.08 milliseconds. The difference came from the desire to compare his model to the code represented in Chapter 7 for the Hodgkin-Huxley model, as well as concerns over the total runtime of the model. Finally, though it is not shown by the figures, the model that we used does not always reach 30 mV. As Izhikevich discusses in his paper, the original graphs were altered to automatically jump to 30 mV when they began to spike to show a uniform graph. This alteration does not show the accuracy of his model very well, but we followed his example and made the necessary changes to the code to guarantee consistency in our results.

Once we had created Figure 2, it was a simple matter to alter the inputted parameters to fit the seven other situations, then use GraphicsGrid to show all eight types of spiking at once in Figure 3.

Figure 3: By repeating the process of creating Figure 2 eight times for a variety of different parameter values, we are able to replicate the model from Izhikevich in the same order as in his paper.

While the previous discrepancies are once again present in this figure, we can say with confidence that the difference between each situation is as prevalent as in the original figure. The first parameter, the time scale of the membrane potential, is closely related to the second parameter, the sensitivity of the recovery variable. As both parameters increase, the ability of the system to return to resting state and spike again increases and allows for more spikes in the same time period as the regularly spiking neuron (such as in the fast spiking and low threshold spiking graphs). When the third parameter, the reset value of the membrane potential, is increased, the neuron also fires faster, but when the effect of the reset value for the recovery variable is simultaneously decreased, the neuron becomes more likely to fire in clusters (such as in the chattering neuron). These results are consistent with the conclusions that Izhikevich came to in his paper, so we can reasonably say that we have replicated his model. Below, we used a Manipulate function to fully illustrate the effect of each parameter on the behavior of the neuron, and invite the reader to follow the link to our code to see the effects firsthand.

Figure 4: In order to allow the reader to explore the conclusions reached about the parameters in this paper, a Manipulate function was used to create another graph. The initial conditions are that of an Regular Spiking Neuron.

Hodgkin-Huxley Model

Having created an accurate representation of Izhikevich's model for a spiking neuron, we now need something to compare it to. We turned to Chapter 7, where the Hodgkin-Huxley model was represented, to find a set of equations and graphs that could be reliably compared to the results we found above. In Figure 5, we have recreated the model using the exact code from Chapter 7, with the exception of an alteration to the injected current and timescale to fit the model from Izhikevich. We did not model all eight different situations because only one is necessary to properly perform a bifurcation analysis.

Figure 5: The Hodgkin-Huxley model shown above has significantly more information than in the paper's model. This was directly made in Chapter 7 of the text, and will be expanded upon in the Bifurcation Analysis.

Basic Comparison of the Two Models

Clearly, the two graphs we have for individual neuron firing in a Regularly Spiking cortical neuron scenario are similar. However, the Hodgkin-Huxley graph has several small dips and hills leading to and coming from the spike that are not present in the very smooth model from the paper. This can be explained by the complexity of their equations; the Izhikevich equations are simplistic in nature because he was attempting to make the model computationally efficient, while the Hodgkin-Huxley equations were designed to be as accurate as possible and had many more parameters as a result. On the surface, it would appear that our hypothesis has been proved. The graphs are remarkably similar, and no significant information appears to have been lost. However, the dips and hills previously observed are concerning, and so we shall also perform a bifurcation analysis on the models to test our hypothesis further.

Extension: Bifurcation Analysis

In order to perform our bifurcation analysis, we need to manipulate our equations into something more reasonable. The way we accomplished this was by referring to the paper again, which cited a different paper by the same author that had a usable equation. The first thing we did was define two sets of parameters, one with a large current added and the other with no current added. Then, we solved the two equations from the second paper using Solve and created two basic nullcline equations to plot. While there are two nullclines, there will be three lines to plot because of our two sets of parameters (the u nullcline will not change from the change in current because it is not dependent on current). When we plot these, we will also plot the VectorPlot of the equations we used, which gives the vector field of the system.

Figure 6: The bifurcation diagram for the paper model. Notice the two equilibrium points in the system when the V nullcline does not have an injected current.

There are several things to note from this system. First, there are two equilibrium points in the 0 current system, which occur at (-60, 0) and (-35, -35). The first is a stable equilibrium point, as the values all head to that value, and the second is unstable. The 70 current system, however, never intersects the u nullcline, which suggests that there are no equilibrium points for that system. If either system is greater than the value of the second, unstable equilibrium point, it will head out to 30 mV before being sent back to its starting point. That is the behavior of the system reaching a spike and then resetting for the next one.

We also needed the bifurcation diagram for the Hodgkin-Huxley model. The following graph comes from a paper which performed this analysis [7].

Figure 7: The bifurcation diagram from Hodgkin-Huxley. The level of detail is much higher than in the paper model. This graph comes from a cited source, as noted above.


Model vs. Hypothesis

In our bifurcation analysis, outlined in the Results section, we were able to demonstrate the behavior of both Izhikevich's original and the Hodgkin-Huxley model demonstrated in Chapter 7 of our textbook. There were considerable differences between the nullclines of the membrane potential, and the value of the equilibrium points of the models were not the same, either. The hypothesis that we were testing was that the Izhikevich did not lose a significant amount of information when he simplified the Hodgkin-Huxley model to be computationally efficient, and, based on the nullclines and equilibrium points, our model does not support our hypothesis. It should be noted that this is in spite of the replications of both models being incredibly visually similar.


There are several inherent limitations of our results that we found difficult or impossible to overcome. First, the two models that we were comparing had a different set of initial parameter values and were constructed differently. The Izhikevich model was meant to be simpler, and so it only had two state variables compared to the four state variables in the Chapter 7 model of the Hodgkin-Huxley equations. These state variables provide the depth of detail that leads to our dismissal of our hypothesis, and yet it is hard to say whether or not the systems would be comparable if we had the same number of state variables.

In addition, as outlined in our model description, there are several assumptions we had to make about both systems that undermine the applicability of our data to real life systems.

Model Discrepancies

Our results section outlined in greater detail the discrepancies between our replications and those originally presented by Izhikevich. In his paper, the code he used artificially spiked the neurons to 30 mV every time to show consistency, while we simply increased the number of data points to have our graph show 30 mV spikes every time. In addition, we used a single set of equations with changing parameters based on what he provided in the paper. However, Izhikevich changed his equations every time he made a graph in order to increase the accuracy of his work. He justified this by pointing out that there are a variety of different neurons out there, and that most instances of modeling do not require being able to address every single one [8]

Literature Review of Conclusions

Once we had reached the conclusion that our hypothesis was not supported by our models, we needed to determine why the models did not match each other. In one of our cited sources from the introduction, the processes that control the neuron's membrane potential include a series of equations governing the exchange of ions across the membrane during firing [1] [2]. The concentration of potassium ions, sodium ions, and other materials changes constantly through the spiking, creating the differences in charge that force the spike. However, in the model created by Izhikevich, these are represented as a series of parameters that are unique to a firing period [4]. Their static state creates graphs that look very similar for each individual scenario, but do not provide the same level of detail as the Hodgkin-Huxley model from Chapter 7, which had a state variable for each of the main ion contributors in the system [5]. These concentration gradients are incredibly important to the membrane potential, and are the most likely reason, according to literature, that the two models have such starkly different nullclines in our bifurcation analysis.

Future Work

A reasonable extension of this work beyond the bifurcation analysis that we conducted would be to see exactly what information was lost between the two models. From our nullcline results, it would appear that the primary loss of detail came in the current that was injected into the system, which could be potentially corrected to provide the biologically plausable and computationally efficient model that Izhikevich wanted to provide. If that proved to be impossible to represent, we could also try to go with our original plan with regards to the comparison, which was to create a network of thousands of spiking cortical neurons in place of a bifurcation analysis to show an increasing level of difference between the systems as the number of neurons increased.


  1. 1.0 1.1 Hagan, Martin; Howard, Demuth; Beale, Mark; De Jess, Orlando (2014). Neural Network Design (2nd ed.). New York: Martin Hagan. 
  2. 2.0 2.1 Rojas, Raul (2013). Neural Networks: A Systematic Introduction (1st ed.). Berlin: Springer. 
  3. John, Staddon; Grossberg, Stephen; Commons, Michael (1991). Quantitative Analyses of Behavior (1st ed.). Hillsdale: Lawrence Elbaum Associates. 
  4. 4.0 4.1 4.2 4.3 4.4 4.5 Izhikevich, Eugene. "Simple Model of Spiking Neurons". IEEE Transactions on Neural Networks. 14 (6): 1569–1572. 
  5. 5.0 5.1 Chiel, Hillel J. (2017). Chapter 7 (18th ed.). Hillel J. Chiel. 
  6. Soleimani, Hamid; Bayandpour, Mohammad; Ahmadi, Arash; Abbott, Derek (April 2017). "Neural Dynamics and Bifurcation". IEEE. 
  7. Guckenheimer, J (September 1993). "Bifurcation of the Hodgkin and Huxley equations: A new twist". Bulletin of Mathematical Biology. 55 (5). 
  8. Izhikevich, E (September 2004). "Which model to use for cortical spiking neurons?". IEEE. 15 (5).