-
Notifications
You must be signed in to change notification settings - Fork 0
/
SA_Francisco.py
64 lines (57 loc) · 2.08 KB
/
SA_Francisco.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import numpy as np
import scipy as sp
import scipy.integrate as sp_int
import matplotlib.pyplot as plt
import matplotlib
# define function that is time derivative of voltage: f = dV/dt
def f(t, y):
# print (t,y, "kek")
# set the relevant constants
C_m, g_K, g_Na, g_l, V_K, V_Na, V_l = x
# set external current: make sure it is 0 for large times.
if t <= 0.00001:
I_e = 1
else:
I_e = 0
# set the variables that are to be integrated
V, n, m, h = y
# define DV/dt ('_dot' denotes time differentiation)
V_dot = 1/C_m * (I_e - ((g_K * (n**4) * (V - V_K) + g_Na * (m**3)*h) *
(V - V_Na) + g_l * (V - V_l)))
# equations governing opening/closing rates.
a_n = 0.01 * (V + 10) / ( np.exp((V+10)/10) - 1)
b_n = 0.125 * np.exp(V/80)
a_m = 0.1 * (V + 25) / (np.exp((V + 25)/10) - 1)
b_m = 4 * np.exp(V / 18)
a_h = 0.07 * np.exp(V / 20)
b_h = 1 / (np.exp((V + 30)/10) + 1)
# enter the equations controlling the gating variables.
# 'g_1' is dn/dt, 'g_2' is dm/dt, 'g_3' is dh/dt
g_1 = a_n * (1 - n) - b_n*n
g_2 = a_m * (1 - m) - b_m*m
g_3 = a_h * (1 - h) - b_h*h
# print("kek2")
# since w = [V, n, m, h] we return [V_dot, n_dot, m_dot, h_dot]
return [V_dot, g_1, g_2, g_3]
# print("hello pig")
# enter the values of the constants. Values taken from Table 3 in `Membrane Current In Nerve`
# values are entered like: conts = [I_e, C_m, g_K, g_Na, g_l, V_K, V_Na, V_l]
conts = [10**(-3), 10**(-3), 10**(-3), 10**(-3), 10**(-3), 10**(-3), -10**(-3)]
# enter intial values for V, n, m, h
V_0 = 10**(-9)
n_0 = 0
m_0 = 0
h_0 = 0
y_0 = [V_0, n_0, m_0, h_0]
# print("pig initialised")
# create timescale. t_interval is the time interval in which to calculate the solution.
# t_points are the points at which the solution is stored.
t_interval = (0.0, 10.0)
numpoints = 1000
t_points = np.linspace(t_interval[0], t_interval[1], numpoints)
x = conts
# solve coupled ODEs with scipy's solver
# print("about to solve pig")
soln = sp_int.solve_ivp(f, t_interval, y_0, 'RK45', t_points)
plt.plot(soln.t, soln.y[0, :])
plt.show()