# set up a dictionary of parameters
E_params = {
'E_leak' : -7.0e-2,
'G_leak' : 3.0e-09,
'C_m' : 3.0e-11,
'I_ext' : 0*1.0e-10
}
Na_params = {
'Na_E' : 5.0e-2,
'Na_G' : 1.0e-6,
'k_Na_act' : 3.0e+0,
'A_alpha_m_act' : 2.0e+5,
'B_alpha_m_act' : -4.0e-2,
'C_alpha_m_act' : 1.0e-3,
'A_beta_m_act' : 6.0e+4,
'B_beta_m_act' : -4.9e-2,
'C_beta_m_act' : 2.0e-2,
'l_Na_inact' : 1.0e+0,
'A_alpha_m_inact' : 8.0e+4,
'B_alpha_m_inact' : -4.0e-2,
'C_alpha_m_inact' : 1.0e-3,
'A_beta_m_inact' : 4.0e+2,
'B_beta_m_inact' : -3.6e-2,
'C_beta_m_inact' : 2.0e-3
}
K_params = {
'k_E' : -9.0e-2,
'k_G' : 2.0e-7,
'k_K' : 4.0e+0,
'A_alpha_m_act' : 2.0e+4,
'B_alpha_m_act' : -3.1e-2,
'C_alpha_m_act' : 8.0e-4,
'A_beta_m_act' : 5.0e+3,
'B_beta_m_act' : -2.8e-2,
'C_beta_m_act' : 4.0e-4
}
params = {
'E_params' : E_params,
'Na_params' : Na_params,
'K_params' : K_params
}Simulating Dynamical Systems II
Modelling Action Potentials
In this chapter we will use a model of voltage-gated ion channels in a single neuron to simulate action potentials. The model is based on the work by Hodgkin & Huxley in the 1940s and 1950s. A good reference to refresh your memory about how ion channels in a neuron work is the Kandel, Schwartz & Jessel book “Principles of Neural Science”.
A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. (Lond.), 117(4):500-544, Aug 1952
A. L. Hodgkin, A. F. Huxley, A. L. Hodgkin, and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. 1952. Bull. Math. Biol., 52(1-2):25-71, 1990
O. Ekeberg, P. Wallen, A. Lansner, H. Traven, L. Brodin, and S. Grillner. A computer based model for realistic simulations of neural networks. I. The single neuron and synaptic interaction. Biol Cybern, 65(2):81-90, 1991
E.R. Kandel, J.H. Schwartz, T.M. Jessell, et al. Principles of neural science, volume 4. McGraw-Hill New York, 2000
Introduction
To model the action potential we will use an article by Ekeberg et al. (1991) published in Biological Cybernetics (see citation above). When reading the article you can focus on the first three pages (up to paragraph 2.3) and try to find answers to the following questions:
How many differential equations are there?
What is the order of the system described in equations 1-9?
What are the states and state derivatives of the system?
Before we begin coding up the model, it may be useful to remind you of a fundamental law of electricity, one that relates electrical potential V to electric current I and resistance R (or conductance G, the reciprocal of resistance). This of course is known as Ohm’s law:
V = IR \tag{1}
or
V = \frac{I}{G} \tag{2}
where V is the voltage, I is the current, and R is the resistance.
Our goal here is to code up a dynamical model of the membrane’s electric circuit including two types of ion channels: sodium and potassium channels. We will use this model to better understand the process underlying the origin of an action potential.
The Neuron Model
The Figure below, adapted from Ekeberg et al., 1991, schematically illustrates the model of a neuron. In panel A we see a soma, and multiple dendrites. Each of these can be modelled by an electrical “compartment” (Panel B) and the passive interactions between them can be modelled as a pretty standard electrical circuit (see Biological Neuron Model for more details about compartmental models of neurons). In panel C, we see an expanded model of the Soma from panel A. Here, a number of active ion channels are included in the model of the soma.
For our purposes here, we will focus on the soma, and we will not include any additional dendrites in our implementation of the model. Thus essentially we will be modelling what appears in panel C, and at that, only a subset.
In panel C we see that the soma can be modelled as an electrical circuit with a sodium ion channel (Na), a potassium ion channel (K), a calcium ion channel (Ca), and a calcium-dependent potassium channel (K(Ca)). What we will be concerned with simulating, ultimately, is the intracellular potential E.
Passive Properties
Equation (1) of Ekeberg is a differential equation describing the relation between the time derivative of the membrane potential E as a function of the passive leak current through the membrane, and the current through the ion channels. Note that Ekeberg uses E instead of the typical V symbol to represent electrical potential.
\frac{dE}{dt} = \frac{(E_{leak}-E)G_{m} + \sum{\left(E_{comp}-E\right)}G_{core} + I_{channels}}{C_{m}} \tag{3}
Don’t panic, it’s not actually that complicated. What this equation is saying is that the rate of change of electrical potential across the current (the left hand side of the equation, \frac{dE}{dt}) is equal to the sum of a bunch of other terms, divided by membrane capacitance C_{m} (the right hand side of the equation). Recall from basic physics that capacitance is a measure of the ability of something to store an electrical charge.
The “bunch of other things” is a sum of three things, actually, (from left to right): a passive leakage current, plus a term characterizing the electrical coupling of different compartments, plus the currents of the various ion channels. Since we are not going to be modelling dendrites here, we can ignore the middle term on the right hand side of the equation \sum{\left(E_{comp}-E\right)}G_{core} which represents the sum of currents from adjacent compartments (we have none).
We are also going to include in our model an external current I_{ext}. This can essentially represent the sum of currents coming in from the dendrites (which we are not explicitly modelling). It can also represent external current injected in a patch-clamp experiment. This is what we as experimenters can manipulate, for example, to see how neuron spiking behaviour changes. So what we will actually be working with is this:
\frac{dE}{dt} = \frac{(E_{leak}-E)G_{m} + I_{channels} + I_{ext}}{C_{m}} \tag{4}
What we need to do now is unpack the I_{channels} term representing the currents from all of the ion channels in the model. Initially we will only be including two, the potassium channel (K) and the sodium channel (Na).
Sodium Channels (Na)
The current through sodium channels that enter the soma are represented by equation (2) in Ekeberg et al. (1991):
I_{Na} = (E_{Na} - E_{soma})G_{Na}m^{3}h \tag{5}
where m is the activation of the sodium channel and h is the inactivation of the sodium channel, and the other terms are constant parameters: E_{Na} is the reversal potential, G_{Na} is the maximum sodium conductance throught the membrane, and E_{soma} is the membrane potential of the soma.
The activation m of the sodium channels is described by the differential equation (3) in Ekeberg et al. (1991):
\frac{dm}{dt} = \alpha_{m}(1-m) - \beta_{m}m \tag{6}
where \alpha_{m} represents the rate at which the channel switches from a closed to an open state, and \beta_{m} is rate for the reverse. These two parameters \alpha and \beta depend on the membrane potential in the soma. In other words the sodium channel is voltage-gated. Equation (4) in Ekeberg et al. (1991) gives these relationships:
\begin{aligned} \alpha_{m} &= &\frac{A(E_{soma}-B)}{1-e^{(B-E_{soma})/C}}\\ \beta_{m} &= &\frac{A(B-E_{soma})}{1-e^{(E_{soma}-B)/C}}\end{aligned} \tag{7}
A tricky bit in the Ekeberg et al. (1991) paper is that the A, B and C parameters above are different for \alpha and \beta even though there is no difference in the symbols used in the equations.
The inactivation of the sodium channels is described by a similar set of equations: a differential equation giving the rate of change of the sodium channel deactivation, from Ekeberg et al. (1991) equation (5):
\frac{dh}{dt} = \alpha_{h}(1-h) - \beta_{h}h \tag{8}
and equations specifying how \alpha_{h} and \beta_{h} are voltage-dependent, given in Ekeberg et al. (1991) equation (6):
\begin{aligned} \alpha_{h} &= &\frac{A(B-E_{soma})}{1-e^{(E_{soma}-B)/C}}\\ \beta_{h} &= &\frac{A}{1-e^{(B-E_{soma})/C}}\end{aligned} \tag{9}
Note again that although the terms A, B and C are different for \alpha_{h} and \beta_{h} even though they are represented by the same symbols in the equations.
So in summary, for the sodium channels, we have two state variables: (m,h) representing the activation (m) and deactivation (h) of the sodium channels. We have a differential equation for each, describing how the rate of change (the first derivative) of these states can be calculated: Ekeberg equations (3) and (5). Those differential equations involve parameters (\alpha,\beta), one set for m and a second set for h. Those (\alpha,\beta) parameters are computed from Ekeberg equations (4) (for m) and (6) (for h). Those equations involve parameters (A,B,C) that have parameter values specific to \alpha and \beta and m and h (see Table 1 of Ekeberg et al., 1991).
Potassium Channels (K)
The potassium channels are represented in a similar way, although in this case there is only channel activation, and no inactivation. In Ekeberg et al. (1991) the three equations (7), (8) and (9) represent the potassium channels:
I_{k} = (E_{k}-E_{soma})G_{k}n^{4} \tag{10}
\frac{dn}{dt} = \alpha_{n}(1-n) - \beta_{n}n \tag{11}
where n is the state variable representing the activation of potassium channels. As before we have expressions for (\alpha,\beta) which represent the fact that the potassium channel is also voltage-gated:
\begin{aligned} \alpha_{n} &= &\frac{A(E_{soma}-B)}{1-e^{(B-E_{soma})/C}}\\ \beta_{n} &= &\frac{A(B-E_{soma})}{1-e^{(E_{soma}-B)/C}}\end{aligned} \tag{12}
Again, the parameter values for (A,B,C) can be found in Ekeberg et al., (1991) Table 1.
To summarize, the potassium channel has a single state variable n representing the activation of the potassium channel.
Summary
We have a model now that includes four state variables:
E representing the potential in the soma, given by differential equation (1) in Ekeberg et al., (1991)
m representing the activation of sodium channels, Ekeberg equation (3)
h representing the inactivation of sodium channels, Ekeberg equation (5)
n representing the activation of potassium channels, Ekeberg equation (8)
Each of the differential equations that define how to compute state derivatives, involve (\alpha,\beta) terms that are given by Ekeberg equations (4) (for m), (6) (for h) and (9) (for n).
So what we have to do in order to simulate the dynamic behaviour of this neuron over time, is simply to implement these equations in MATLAB code, give the system some reasonable initial conditions, and simulate it over time using the ode45() function.
Python code
We begin by setting up all of the necessary model parameters (there are many). They are found in Ekeberg et al. (1991) Tables 1 and 2. I have chosen to do this using a Python dictionary. This is a useful data type to parcel all of our parameters together. Unlike an array or list, which we would have to index using integer values (and then keep track of which one corresponded to which parameter), with a dictionary, we can index into it using string labels.
Our next bit of code is the implementation of the ODE function itself. Remember, our ultimate goal is to model E, the potential across the soma membrane. We know from Ekeberg equation (1) that the rate of change of E depends on leakage current and on the sum of currents from other channels. These other currents are given by Ekeberg equations (2) (sodium) and (7) (potassium). These equations involve the three other states in our system: sodium activation m, sodium inactivation h and potassium activation n, which are each defined by their own differential equations, Ekeberg equations (3), (5) and (8), respectively. It’s just a matter of coding things up step by step.
# define our ODE function
def neuron(state, t, params):
"""
Purpose: simulate Hodgkin and Huxley model for the action potential using
the equations from Ekeberg et al, Biol Cyb, 1991.
Input: state ([E m h n] (ie [membrane potential; activation of
Na++ channel; inactivation of Na++ channel; activation of K+
channel]),
t (time),
and the params (parameters of neuron; see Ekeberg et al).
Output: statep (state derivatives).
"""
E = state[0]
m = state[1]
h = state[2]
n = state[3]
Epar = params['E_params']
Na = params['Na_params']
K = params['K_params']
# external current (from "voltage clamp", other compartments, other neurons, etc)
I_ext = Epar['I_ext']
# calculate Na rate functions and I_Na
alpha_act = Na['A_alpha_m_act'] * (E-Na['B_alpha_m_act']) / (1.0 - np.exp((Na['B_alpha_m_act']-E) / Na['C_alpha_m_act']))
beta_act = Na['A_beta_m_act'] * (Na['B_beta_m_act']-E) / (1.0 - np.exp((E-Na['B_beta_m_act']) / Na['C_beta_m_act']) )
dmdt = ( alpha_act * (1.0 - m) ) - ( beta_act * m )
alpha_inact = Na['A_alpha_m_inact'] * (Na['B_alpha_m_inact']-E) / (1.0 - np.exp((E-Na['B_alpha_m_inact']) / Na['C_alpha_m_inact']))
beta_inact = Na['A_beta_m_inact'] / (1.0 + (np.exp((Na['B_beta_m_inact']-E) / Na['C_beta_m_inact'])))
dhdt = ( alpha_inact*(1.0 - h) ) - ( beta_inact*h )
# Na-current:
I_Na =(Na['Na_E']-E) * Na['Na_G'] * (m**Na['k_Na_act']) * h
# calculate K rate functions and I_K
alpha_kal = K['A_alpha_m_act'] * (E-K['B_alpha_m_act']) / (1.0 - np.exp((K['B_alpha_m_act']-E) / K['C_alpha_m_act']))
beta_kal = K['A_beta_m_act'] * (K['B_beta_m_act']-E) / (1.0 - np.exp((E-K['B_beta_m_act']) / K['C_beta_m_act']))
dndt = ( alpha_kal*(1.0 - n) ) - ( beta_kal*n )
I_K = (K['k_E']-E) * K['k_G'] * n**K['k_K']
# leak current
I_leak = (Epar['E_leak']-E) * Epar['G_leak']
# calculate derivative of E
dEdt = (I_leak + I_K + I_Na + I_ext) / Epar['C_m']
statep = [dEdt, dmdt, dhdt, dndt]
return statepNext we run a simulation by setting up our initial states, and a time array, and then calling odeint(). Note that we are injecting some external current by changing the value of the params['E_params']['I_ext'] entry in the params dictionary.
import numpy as np
import scipy as sp
# simulate
# set initial states and time vector
state0 = np.array([-70e-03, 0, 1, 0])
t = np.arange(0, 0.2, 0.001)
# let's inject some external current
params['E_params']['I_ext'] = 1.0e-10
# run simulation
state = sp.integrate.odeint(neuron, state0, t, args=(params,))Finally, we plot the results:
import matplotlib.pyplot as plt
plt.figure(figsize=(6,9))
plt.subplot(4,1,1)
plt.plot(t, state[:,0])
plt.title('membrane potential')
plt.subplot(4,1,2)
plt.plot(t, state[:,1])
plt.title('Na2+ channel activation')
plt.subplot(4,1,3)
plt.plot(t, state[:,2])
plt.title('Na2+ channel inactivation')
plt.subplot(4,1,4)
plt.plot(t, state[:,3])
plt.title('K+ channel activation')
plt.xlabel('TIME (sec)')
plt.tight_layout()
plt.show()
Things To Try
Alter the
neuron()function so that the modelled neuron only has the leakage current and external current. In other words, comment out the terms related to sodium and potassium channels. Run a simulation with an initial membrane potential of -70 mv and an external current of 0.0 mv. What happens and why?Change the external current to 1.0e-10 and re-run the simulation. What happens and why?
Add in the terms related to the sodium channel (activation and deactivation). Run a simulation with external current of 1.0e-10 and initial states [-70e-03, 0, 1]. What happens and why?
Add in the terms related to the potassium channel. Run a simulation with external current of 1.0e-10 and initial states [-70e-03, 0, 1, 0]. What happens and why?
Play with the external current level (increase it slightly, decrease it slightly, etc). What is the effect on the behaviour of the neuron?
What is the minimum amount of external current necessary to generate an action potential? Why?
MATLAB version
Here is a link to a MATLAB version of these notes from a previous version of the course.
thanks to Dinant Kistemaker for co-authoring this topic