<jupyter_start><jupyter_text>--- Neuron modelA *membrane equation* and a *reset condition* define our *leaky-integrate-and-fire (LIF)* neuron:\begin{align*}\\&\tau_m\,\frac{d}{dt}\,V(t) = E_{L} - V(t) + R\,I(t) &\text{if }\quad V(t) \leq V_{th}\\\\&V(t) = V_{reset} &\text{otherwise}\\\\\end{align*}where $V(t)$ is the membrane potential, $\tau_m$ is the membrane time constant, $E_{L}$ is the leak potential, $R$ is the membrane resistance, $I(t)$ is the synaptic input current, $V_{th}$ is the firing threshold, and $V_{reset}$ is the reset voltage. We can also write $V_m$ for membrane potential - very convenient for plot labels.The membrane equation is an *ordinary differential equation (ODE)* that describes the time evolution of membrane potential $V(t)$ in response to synaptic input and leaking of change across the cell membrane.**Note that, in this tutorial the neuron model will not implement a spiking mechanism.**<jupyter_code>#-------------------------------Leaky Integrate-and-Fire (LIF) neuron----------------------------------------


# Import libraries
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import YouTubeVideo


# defining and initializing the main simulation variables
t_max = 150e-3   # second
dt = 1e-3        # second
tau = 20e-3      # second
el = -60e-3      # milivolt
vr = -70e-3      # milivolt
vth = -50e-3     # milivolt
r = 100e6        # ohm
i_mean = 25e-11  # ampere

print(t_max, dt, tau, el, vr, vth, r, i_mean)
#-----------------------------simulate the synaptic input I(t) with a sinusoidal model---------------------------------
# initialize t
t = 0

# loop for 10 steps, variable 'step' takes values from 0 to 9
for step in range(10):
  t = step * dt
  # --> insert your code here
  i = i_mean * (1 + np.sin(2 * np.pi * t / 0.01)) 
  print(f'{t:.3f} {i:.4e}')
#-----------------------------ODE integration without spikes------------------------------------------------------------
# initialize step_end and v
step_end = 10
v = el

# loop for step_end steps

# --> insert your code here
for step in range(step_end):
  t = step * dt
  i = i_mean * (1 + np.sin(2 * np.pi * t / 0.01)) 
  print(f"{t:.3f} {v:.4e}")      # print before update!!!!
  v = v + (el- v + r * i) * dt / tau
#---------------------------Plot the values of I(t) between t=0 and t=0.024------------------
# initialize step_end
step_end = 25

# initialize the figure
plt.figure()

# loop for step_end steps
for step in range(step_end):
  t = step * dt
  i = i_mean * (1 + np.sin(2 * np.pi * t / 0.01)) 
  plt.plot(t,i,'ko')

plt.title('Synaptic Input $I(t)$')
plt.xlabel('time (s)')
plt.ylabel(r'$I$ (A)')
plt.show()
#------------------------Plot the values of V(t) between t=0 and t=tmax---------------------------
# initialize step_end and v
step_end = int(t_max/dt)   # t_max = step * dt
v = el

for step in range(step_end):
  t = step * dt
  i = i_mean * (1 + np.sin(2 * np.pi * t / 0.01)) 
  plt.plot(t,v,'k.')  # plot before update!!!!
  v = v + (el- v + r * i) * dt / tau

t = t + dt
plt.plot(t, v, 'k.')

plt.title('$V_m$ with sinusoidal I(t)')
plt.xlabel('time (s)')
plt.ylabel(r'$V_m$ (V)')
plt.show()
#-------------------------Plot the values of V(t) between t=0 and t=tmax−Δt with random input I(t)--------------------------
# set random number generator
np.random.seed(2020)

# initialize step_end and v
step_end = int(t_max / dt)
v = el

# initialize the figure
plt.figure()
plt.title('$V_m$ with random I(t)')
plt.xlabel('time (s)')
plt.ylabel(r'$V_m$ (V)')

# loop for step_end steps
for step in range(step_end):
  t = step * dt
  i = i_mean * (1 + 0.1 * (t_max/dt)**0.5 * (2*np.random.random()-1))
  plt.plot(t,v,'k.')
  v = v + (el- v + r * i) * dt / tau

plt.show()
#------------------------Plot multiple realizations (N=50) of V(t) by storing in a list the voltage of each neuron at time t--------------------
# set random number generator
np.random.seed(2020)

# initialize step_end, n and v_n
step_end = int(t_max / dt)
n = 50
v_n = [el] * n   # python中 列表*数字 就是将列表里的数字重复n遍，存储到新的列表中,不改变原来列表里的值

# initialize the figure
plt.figure()
plt.title('Multiple realizations of $V_m$')
plt.xlabel('time (s)')
plt.ylabel('$V_m$ (V)')

# loop for step_end steps
for step in range(step_end):
  t = step * dt
  plt.plot([t]*n,v_n,'k.', alpha=0.05)

  for j in range(n):
    i = i_mean * (1 + 0.1 * (t_max/dt)**0.5 * (2*np.random.random()-1))
    v_n[j] = v_n[j] + (el- v_n[j] + r * i) * dt / tau

plt.show()
#------------------------Add the sample mean to the plot-------------------------------------------------------------
# set random number generator
np.random.seed(2020)

# initialize step_end, n and v_n
step_end = int(t_max / dt)
n = 50
v_n = [el] * n

# initialize the figure
plt.figure()
plt.title('Multiple realizations of $V_m$')
plt.xlabel('time (s)')
plt.ylabel('$V_m$ (V)')

# loop for step_end steps
for step in range(step_end):
  t = step * dt
  v_mean = sum(v_n)/n
  plt.plot([t] * n, v_n, 'k.', alpha=0.05)
  plt.plot(t,v_mean, 'C0.', alpha=0.8, markersize=10)

  # loop for n steps
  for j in range(0, n):
    i = i_mean * (1 + 0.1 * (t_max/dt)**(0.5) * (2* np.random.random() - 1))
    #if j == 0 :
     # continue   这里按理是要加上这一步的
    v_n[j] = v_n[j] + (dt / tau) * (el - v_n[j] + r*i)
    
 
plt.show()
#----------------------------Add the sample standard deviation to the plot---------------------------------------
# set random number generator
np.random.seed(2020)

# initialize n, v_n and step_end
step_end = int(t_max / dt)
n = 50
v_n = [el] * n

# initialize the figure
plt.figure()
plt.title('Multiple realizations of $V_m$')
plt.xlabel('time (s)')
plt.ylabel('$V_m$ (V)')

# loop for step_end steps
for step in range(step_end):
  t = step * dt
  v_mean = sum(v_n)/n
  v_var_n = [(v - v_mean)**2 for v in v_n]  #不能直接写成 v_n - v_mean，也就是列表和数字不能直接相减, 但是如果用numpy.array的数字就是elementwise操作，就可以直接加减乘除
  v_var = sum(v_var_n) / n
  v_std = np.sqrt(v_var)

  plt.plot([t] * n, v_n, 'k.', alpha=0.05)
  plt.plot(t,v_mean, 'C0.', alpha=0.8, markersize=10)
  plt.plot(t,v_mean + v_std,'C7.', alpha=0.8)
  plt.plot(t,v_mean - v_std,'C7.', alpha=0.8)

  # loop for n steps
  for j in range(0, n):
    i = i_mean * (1 + 0.1 * (t_max/dt)**(0.5) * (2* np.random.random() - 1))
    v_n[j] = v_n[j] + (dt / tau) * (el - v_n[j] + r*i)
#------------------------------Rewrite the single neuron plot with random input from Exercise 7 with numpy arrays. The time range, voltage values, and synaptic current are initialized or pre-computed as numpy arrays before numerical integration.
# Exercise 11
# set random number generator
np.random.seed(2020)

# initialize step_end, t_range, v and syn
step_end = int(t_max/dt)
t_range = np.linspace(0, t_max, num=step_end)  # linspce在指定的间隔内返回均匀间隔的数字，而之前的写法是range(step_end)控制循环次数，在每次循环里一个点一个点的打印t，并更新t=t*dt
                                               # 这里说明不需要一个点一个点的去打印了，而是将t全部存储为np.array
v = el * np.ones(step_end)    #对比列表: v_n = [el]*n
syn = i_mean * (1 + 0.1 * (t_max/dt) ** (0.5) * (2 * np.random.random(step_end) - 1))  #现在i也是一次性存np.array

# initialize the figure
plt.figure()
plt.title('$V_m$ with random I(t)')
plt.xlabel('time (s)')
plt.ylabel(r'$V_m$ (V)')

# loop for step_end steps
#for step in range(step_end):
  #if step == 0:
   # continue
  #v[step] = v[step] + (el- v[step] + r * syn[step]) * dt / tau    #这样写是错误的！！！因为这里一开始直接将每个时间点t上的v都初始化为el，
                                                                   # 而在前面由于是一个一个打印的，所以当前的t就是上一次的t，每一次都将t给覆盖掉，现在我们用np.array去存储下来

# loop for step_end - 1 steps
for step in range(1, step_end):
  v[step] = v[step - 1] + (dt / tau) * (el - v[step - 1] + r * syn[step])

plt.plot(t_range,v,'k.')
plt.show()
#-----------------------------Let's practice using enumerate to iterate over the indexes and values of the synaptic current array syn----------------------------
#普通的for循环 
#i = 0 
#seq = ['one', 'two', 'three'] 
#for element in seq:     
#    print i, seq[i]     
#    i +=1 

#enumerate( ) 
#seq = ['one', 'two', 'three'] 
#for i, element in enumerate(seq):      
#    print i, element

# set random number generator
np.random.seed(2020)

# initialize step_end, t_range, v and syn
step_end = int(t_max / dt)
t_range = np.linspace(0, t_max, num=step_end)
v = el * np.ones(step_end)
syn = i_mean * (1 + 0.1 * (t_max / dt)**(0.5) * (2 * np.random.random(step_end) - 1))

# loop for step_end values of syn
for step, i in enumerate(syn):

  # skip first iteration
  if step==0:
    continue

  v[step] = v[step - 1] + (dt / tau) * (el - v[step - 1] + r*i)

with plt.xkcd():
  # initialize the figure
  plt.figure()
  plt.title('$V_m$ with random I(t)')
  plt.xlabel('time (s)')
  plt.ylabel(r'$V_m$ (V)')
  
  plt.plot(t_range, v, 'k')
  plt.show()
#-------------------------Plot multiple realizations (N=50) of V(t) by storing the voltage of each neuron at time t in a numpy array---------------------------
# set random number generator
np.random.seed(2020)

# initialize step_end, n, t_range, v and syn
step_end = int(t_max / dt)
n = 50
t_range = np.linspace(0, t_max, num=step_end)
v_n = el * np.ones([n, step_end])     # 存成二维数组
syn = i_mean * (1 + 0.1 * (t_max / dt)**(0.5) * (2 * np.random.random([n, step_end]) - 1))   #二维数组里都是随机数

for step in range(1, step_end):
  v_n[:,step] = v_n[:,step - 1] + (dt / tau) * (el - v_n[:,step - 1] + r * syn[:,step])  #syn即i的不一样，造成对于某个时间点t,50次重复中它的v不一样


# initialize the figure
plt.figure()
plt.title('$V_m$ with random I(t)')
plt.xlabel('time (s)')
plt.ylabel(r'$V_m$ (V)')

plt.plot(t_range,v_n.T,'k', alpha=0.3)   #这里必须是v_n.T，因为它是二维数组，相当于行数要与t_range相同，才能映射到列数v

plt.show()
#---------------------------Add sample mean and standard deviation to the plot-----------------------------------------------------------
# set random number generator
np.random.seed(2020)

# initialize step_end, n, t_range, v and syn
step_end = int(t_max / dt)
n = 50
t_range = np.linspace(0, t_max, num=step_end)
v_n = el * np.ones([n, step_end])
syn = i_mean * (1 + 0.1 * (t_max / dt)**(0.5) * (2 * np.random.random([n, step_end]) - 1))

# loop for step_end - 1 steps
for step in range(1, step_end):

  v_n[:,step] = v_n[:,step - 1] + (dt / tau) * (el - v_n[:, step - 1] + r * syn[:, step])

v_mean = np.mean(v_n,axis = 0)
v_std = np.std(v_n,axis = 0)

plt.figure()
plt.title('Multiple realizations of $V_m$')
plt.xlabel('time (s)')
plt.ylabel('$V_m$ (V)')

plt.plot(t_range, v_n[:-1].T, 'k', alpha=0.3)
plt.plot(t_range, v_n[-1], 'k', alpha=0.3, label='V(t)')
plt.plot(t_range, v_mean, 'C0', alpha=0.8, label='mean')
plt.plot(t_range, v_mean+v_std, 'C7', alpha=0.8)
plt.plot(t_range, v_mean-v_std, 'C7', alpha=0.8, label='mean $\pm$ std')

  # # alternative for filling plot between v_mean-+v_std
  # plt.fill_between(t_range, v_mean-v_std, v_mean+v_std, color='C0', alpha=0.6)
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22869830696641688, 'fraction_numerical': 0.02285132759737256, 'mean_word_length': 2.5664797096667766, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '318087', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4572, 'n_tokens_neox': 4377, 'n_words': 1169, 'path': 'LIF_neuron.ipynb', 'repo_name': 'liuxiaomiao123/NeuroMathAcademy', 'size': '10809'}