%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
import numpy as np
import numpy.linalg as la
from scipy.integrate import RK45
import sympy as sym
from table_maker import *
from functools import partial
from itertools import *
from math import ceil
def cos_bell(x, center=0, width=2*np.pi, height=1):
return (np.cos((x-center)/width*2*np.pi)+1)/2*height * np.heaviside(x-center+width/2,0) * np.heaviside(-x+center+width/2,0)
import matplotlib
matplotlib.rcParams.update({'font.size': 24})
From Kilpatrick and Bressloff 2010 $$\begin{align*} \mu u_t &= -u + \int_{-\infty}^\infty w(x,x^\prime) q(x^\prime,t) f( u(x^\prime,t) - a(x^\prime,t)) \ dx^\prime \\ q_t &= \frac{1 - q}{\alpha} - \beta q f(u - a) \\ \epsilon a_t &= -a + \gamma f(u - a) \end{align*}$$
Modified version - remove synaptic depression $q$ $$\begin{align*} \mu u_t &= -u + \int_{-\infty}^\infty w(x,x^\prime) f( u(x^\prime,t) - a(x^\prime,t)) \ dx^\prime \\ \alpha a_t &= -a + \gamma f(u - a) \end{align*}$$ (note that parameters have been relabeled)
###############################################
# Parameters ##################################
θ = .1
α = 5
γ = 1
μ = 1
def firing_rate(u):
return np.heaviside(u-θ, .5)
def w(x,y):
return .5*np.exp( - np.abs( np.subtract.outer(x, y) ) )
###############################################
# Simulation parameters #######################
# space
a, b = -100, 300
n = 10**3 * 1
xs = np.linspace(a,b,n)
# time
t0 = 0
t_final = 50
k = 1e-2
# initial conditions
# test 1
u0 = np.zeros((2,n))
u0[0] = cos_bell(xs, center=0, width=50)
# Time integrator
def RK4_step(F, t, u, dt):
k1 = F(t,u)
k2 = F(t+dt/2, u + dt/2*k1)
k3 = F(t+dt/2, u + dt/2*k2)
k4 = F(t+dt, u + dt*k3)
return u + dt/6*(k1+2*k2+2*k3+k4)
def F(t,u):
temp = firing_rate(u[0] - u[1])
ret = np.array([
1/μ * (-u[0] + M@(temp)),
(-u[1] + γ*temp)/α
])
return ret
###############################################
# Simulation ##################################
h = xs[1]-xs[0]
M = w(xs, xs) * h
########################
steps = int(np.ceil( (t_final - t0)/k ))
k = (t_final - t0)/steps
us = [u0]
ts = [t0]
for step in range(steps):
ts += [ts[-1] + k]
us += [ RK4_step(F, ts[-1], us[-1], k) ]
print('step %d/%d' % (step,steps), end='\r')
js = [ us[i][0]-us[i][1] for i in range(steps+1) ]
# approximate pulse-width and speed
width = ( np.sum(js[-1]>=θ) - 1 )*h # number of intervals above θ
speed_t0, speed_tf = 20, 40
speed_t0_index, speed_tf_index = [int(t/k) for t in [speed_t0, speed_tf]]
pulse_start = np.max(xs[us[speed_t0_index][0]>θ])
pulse_stop = np.max(xs[us[speed_tf_index][0]>θ])
speed = (pulse_stop - pulse_start)/(ts[speed_tf_index] - ts[speed_t0_index])
print('Width: %g' % width)
print('Speed: %g' % speed)
import pickle
params = (width, speed, γ, α, θ)
pickle.dump(params, open('parameters.pickle', 'wb'))
###############################################
# Animation ###################################
stride = 20
fig, axes = plt.subplots(2, 1, figsize=(12,12))
line_u, = axes[0].plot(xs, us[0][0], 'b-', label="$u$")
line_a, = axes[0].plot(xs, us[0][1], 'm-', label="$a$")
line_j, = axes[1].plot(xs, js[0], 'g-', label="$j=u-a$")
axes[1].plot(xs, θ + 0*xs, 'k--', label='$\\theta$')
y_min = np.min(us)
y_max = np.max(us)
y_min -= .05*np.abs(y_max - y_min)
y_max += .05*np.abs(y_max - y_min)
window = y_min, y_max
axes[0].set_ylim(*window)
y_min = np.min(js)
y_max = np.max(js)
y_min -= .05*np.abs(y_max - y_min)
y_max += .05*np.abs(y_max - y_min)
window = y_min, y_max
axes[1].set_ylim(*window)
for ax in axes:
ax.legend(loc='right')
ax.set_xlim(-50, b)
# ax.set_xlim(a,b)
Δ = -np.log(1 - 4*θ)
plt.plot([150, 150+Δ], [θ]*2, 'k--')
def animate(i):
print('step %d/%d' % (i,len(ts)), end='\r')
line_u.set_ydata(us[i][0])
line_a.set_ydata(us[i][1])
line_j.set_ydata(js[i])
return line_u,
# Init only required for blitting to give a clean slate.
def init():
line_u.set_ydata(us[0])
return line_u,
anim = animation.FuncAnimation(fig, animate, np.arange(0,len(ts),stride), init_func=init,
interval=1/24*1000, blit=True)
plt.close()
HTML(anim.to_html5_video())
###############################################
# Parameters ##################################
θ = .1
α = 5
γ = 1
μ = 1
def firing_rate(u):
return np.heaviside(u-θ, .5)
def w(x,y):
return .5*np.exp( - np.abs( np.subtract.outer(x, y) ) )
###############################################
# Simulation parameters #######################
# space
a, b = -100, 600
n = 10**3 * 1
xs = np.linspace(a,b,n)
# time
t0 = 0
t_final = 100
k = 1e-2
# initial conditions
# test 1
u0 = np.zeros((2,n))
u0[0] = cos_bell(xs, center=0, width=50) + cos_bell(xs, center=500, width=50)
# Time integrator
def RK4_step(F, t, u, dt):
k1 = F(t,u)
k2 = F(t+dt/2, u + dt/2*k1)
k3 = F(t+dt/2, u + dt/2*k2)
k4 = F(t+dt, u + dt*k3)
return u + dt/6*(k1+2*k2+2*k3+k4)
def F(t,u):
temp = firing_rate(u[0] - u[1])
ret = np.array([
1/μ * (-u[0] + M@(temp)),
(-u[1] + γ*temp)/α
])
return ret
###############################################
# Simulation ##################################
h = xs[1]-xs[0]
M = w(xs, xs) * h
########################
steps = int(np.ceil( (t_final - t0)/k ))
k = (t_final - t0)/steps
us = [u0]
ts = [t0]
for step in range(steps):
ts += [ts[-1] + k]
us += [ RK4_step(F, ts[-1], us[-1], k) ]
print('step %d/%d' % (step,steps), end='\r')
###############################################
# Animation ###################################
y_min = np.min(us)
y_max = np.max(us)
y_min -= .05*np.abs(y_max - y_min)
y_max += .05*np.abs(y_max - y_min)
window = y_min, y_max
stride = 30
fig, ax = plt.subplots(1, 1, figsize=(12,8))
line_u, = ax.plot(xs, us[0][0], 'b-', label="$u$")
line_a, = ax.plot(xs, us[0][1], 'm-', label="$a$")
ax.legend(loc='right')
ax.set_ylim(*window)
ax.set_xlim(-50, 550)
# ax.set_xlim(a,b)
def animate(i):
print('step %d/%d' % (i,len(ts)), end='\r')
line_u.set_ydata(us[i][0])
line_a.set_ydata(us[i][1])
return line_u,
# Init only required for blitting to give a clean slate.
def init():
line_u.set_ydata(us[0])
return line_u,
anim = animation.FuncAnimation(fig, animate, np.arange(0,len(ts),stride), init_func=init,
interval=1/24*1000, blit=True)
plt.close()
HTML(anim.to_html5_video())