In Shriki et al., the following equations and parameters were used:
d
x
d
t
=
x
∞

x
τ
x

x
=
h
,
n
,
b
d
x
d
t
=
x
∞

x
τ
x

x
=
h
,
n
,
b
(6)
x
∞
=
α
x
α
x

β
x

x
=
h
,
n
,
m
x
∞
=
α
x
α
x

β
x

x
=
h
,
n
,
m
(7)
τ
x
=
φ
α
x
+
β
x

x
=
h
,
n
,
m
.
τ
x
=
φ
α
x
+
β
x

x
=
h
,
n
,
m
.
(8)
α
m
=

0
.
1
V
+
30
exp
(

0
.
1
(
V
+
30
)
)

1
,
β
m
=
4
exp
(

(
V
+
55
)
/
18
)
,
α
m
=

0
.
1
V
+
30
exp
(

0
.
1
(
V
+
30
)
)

1
,
β
m
=
4
exp
(

(
V
+
55
)
/
18
)
,
(9)
α
h
=
0
.
07
exp
(

(
V
+
44
)
/
20
)
,
β
h
=
1
(
exp
(

0
.
1
(
V
+
14
)
)
+
1
,
α
h
=
0
.
07
exp
(

(
V
+
44
)
/
20
)
,
β
h
=
1
(
exp
(

0
.
1
(
V
+
14
)
)
+
1
,
(10)
α
n
=

0
.
01
(
V
+
34
)
(
exp
(

0
.
1
(
V
+
34
)
)

1
,
β
n
=
0
.
125
exp
(

(
V
+
44
)
/
80
)
,
α
n
=

0
.
01
(
V
+
34
)
(
exp
(

0
.
1
(
V
+
34
)
)

1
,
β
n
=
0
.
125
exp
(

(
V
+
44
)
/
80
)
,
(11)
a
∞
=
1
exp
(

(
V
+
50
)
/
20
)
+
1
,
b
∞
=
1
exp
(
(
V
+
80
)
/
6
)
+
1
.
a
∞
=
1
exp
(

(
V
+
50
)
/
20
)
+
1
,
b
∞
=
1
exp
(
(
V
+
80
)
/
6
)
+
1
.
(12)We first solve for the initial condition V(0)V(0) by setting V'(t)=0V'(t)=0 in Equation 1 since it is generally assumed that the cell is in its steady state configuration at time t=0t=0.
Then we solve Equation 6 using Backwards Euler, noting that Equation 9, Equation 10, Equation 11, Equation 12
use VV from the previous time step. The result allows us to solve Equation 1 using Backwards Euler, govern by the equation
V
j
+
1
=
C
m
Δ
t
V
j
+
g
L
E
L
+
g
¯
N
a
m
∞
3
h
E
N
a
+
g
¯
k
n
4
E
k
+
g
¯
A
a
∞
3
b
E
k
C
m
Δ
t
+
g
L
+
g
¯
N
a
m
∞
3
h
+
g
¯
k
n
4
+
g
¯
A
a
∞
3
b
.
V
j
+
1
=
C
m
Δ
t
V
j
+
g
L
E
L
+
g
¯
N
a
m
∞
3
h
E
N
a
+
g
¯
k
n
4
E
k
+
g
¯
A
a
∞
3
b
E
k
C
m
Δ
t
+
g
L
+
g
¯
N
a
m
∞
3
h
+
g
¯
k
n
4
+
g
¯
A
a
∞
3
b
.
(13)In MATLAB, Equation 13 is written as
top = v(i) + dt/Cm*(g_L*E_L+gbar_na*m_inf^3*h(i+1)*E_na+...
gbar_k*n(i+1)^4*E_k+gbar_A*a_inf^3*b(i+1)*E_k+I_app);
bottom = 1 + dt/Cm*(g_L+...
gbar_na*m_inf^3*h(i+1)+gbar_k*n(i+1)^4+gbar_A*a_inf^3*b(i+1));
v(i+1) = top/bottom;
We say that a cell has spiked, or generated an action potential, if V>30mV.V>30mV. By varying the injective current (kept constant during each run) per simulation, we count the number of spikes in that simulation. This produces an fI curve.