# 002: testi in inglese
import numpy as np
#import sys
import matplotlib.pyplot as plt
import os.path
import re
import time
def AGLIFsynaptic(neuronParameters,equilibriumParameters,currentFileName,tSpikeOutputFileName,voltageOutputFileName):
EL = neuronParameters[0]
vres = neuronParameters[1]
vtm = neuronParameters[2]
Cm = neuronParameters[3]
ith = neuronParameters[4]
tao_m = neuronParameters[5]
sc = neuronParameters[6]
bet = neuronParameters[7]
delta1 = neuronParameters[8]
cost_idep_ini = neuronParameters[9]
Idep_ini_vr = neuronParameters[10]
psi1 = neuronParameters[11]
a=neuronParameters[12]
b=neuronParameters[13]
c=neuronParameters[14]
alp=neuronParameters[15]
istim_min_spikinig_exp=neuronParameters[16]
istim_max_spikinig_exp=neuronParameters[17]
time_scale = 1 / (-sc / (Cm * EL))
H = (90+EL)*sc*(bet-delta1)/(EL*(-200))
sim_lenght = neuronParameters[18]
retteParParsed = neuronParameters[19]
# equilibrium params
v_min = equilibriumParameters[0]
minCurr = equilibriumParameters[1]
zeta = equilibriumParameters[2]
eta = equilibriumParameters[3]
rho = equilibriumParameters[4]
csi = equilibriumParameters[5]
sampling = 2
d_dt = 0.05*sampling
corr_list = []
times_list = []
times_list = np.loadtxt(currentFileName, usecols=(0),dtype=np.float64, unpack=True)
corr_list = np.loadtxt(currentFileName, usecols=(1),dtype=np.float64, unpack=True)
cor = np.array(corr_list)
cor = cor[0:len(cor):1]
times = np.array(times_list)
times = times[0:len(times):1]
def tagliorette(corr,retteParParsed):
vinc_sup = retteParParsed[0]
coeffSup = retteParParsed[1]
constSup = retteParParsed[2]
vinc_inf = retteParParsed[3]
coeffInf = retteParParsed[4]
constInf = retteParParsed[5]
dur_sign=np.inf
if corr<vinc_inf and corr>0:
dur_sign = coeffInf*corr + constInf
if corr>vinc_sup:
dur_sign = coeffSup*corr + constSup
return dur_sign
def V(t, delta, Psi, alpha, beta, IaA0, IdA0, t0, V0):
return (1 / 2) * (beta + (-1) * delta) ** (-1) * (beta ** 2 + ((-1) + beta) * delta) ** (-1) * (4 * beta + (- 1) * (1 + delta) ** 2) ** (-1) * Psi * (2 * np.exp(((-1) * t + t0) * beta) * IdA0 * ((-1) + beta) * beta * (beta + (-1) * delta) * Psi + (-2) * (alpha + (-1) * beta + delta) * (beta ** 2 + ((- 1) + beta) * delta) * Psi + np.exp((1 / 2) * (t + (-1) * t0) * ((-1) + delta + (-1) * Psi)) * (IdA0 * beta * (beta + (-1) * delta) * ((-1) + (-1) * delta + beta * (3 + delta + (-1) * Psi) + Psi) + (-1) * (beta ** 2 + (-1) * delta + beta * delta) * (alpha * (1 + (-2) * beta + delta + (-1) * Psi) + (beta + (-1) * delta) * ((-1) + 2 * IaA0 * beta + (-1) * delta + Psi + V0 * ((-1) + (-1) * delta + Psi)))) + np.exp((1 / 2) * (t + (-1) * t0) * ((-1) + delta + Psi)) * ((-1) * IdA0 * beta * (beta+(-1) * delta) * ((-1) + (-1) * delta + (-1) * Psi + beta * (3 + delta + Psi)) + (beta ** 2 + (-1) * delta+beta * delta) * (alpha * (1 + (-2) * beta + delta + Psi) + (beta + (-1) * delta) * ((-1) + 2 * IaA0 * beta+(-1) * delta + (-1) * Psi + (-1) * V0 * (1 + delta + Psi)))))
def Iadap(t, delta, Psi, alpha, beta, IaA0, IdA0, t0, V0):
return (-2*alpha*(-4*beta**3+beta**2*(-1+delta)**2-delta*(1+delta)**2+beta*delta*(5+2*delta+delta**2))+2*np.exp(((-1)*t + t0) * beta)*IdA0*beta*(4*beta**2+delta*(1+delta)**2-beta*(1+6*delta+delta**2))+np.exp((1 / 2)*(t-t0)*(-1+delta+Psi))*(-IdA0*beta*(beta-delta)*(-1+(-2)*delta-delta**2-Psi+delta*Psi+2*beta*(2+Psi))+(beta**2-delta+beta*delta)*(alpha*(1+(-4)*beta+2*delta+delta**2+Psi-delta*Psi)+(beta-delta)*(4*IaA0*beta-2*(1+V0)*Psi+IaA0*(1+delta)*(-1-delta+Psi))))+np.exp((-1)*(1 / 2) * (t-t0) * (1-delta+Psi))*(IdA0*beta*(beta-delta)*(1+2*delta+delta**2-Psi+delta*Psi+2*beta*(-2+Psi))+(beta**2-delta+beta*delta)*(alpha*(1-4*beta+2*delta+delta**2-Psi+delta*Psi)-(beta-delta)*(-4*IaA0*beta-2*(1+V0)*Psi+IaA0*(1+delta)*(1+delta+Psi)))))/(2*(beta-delta)*(beta**2+(-1+beta)*delta)*(4*beta-(1+delta)**2))
def Idep(t, beta, IdA0, t0):
return np.exp(((-1) * t + t0) * beta) * IdA0
def exp_cum(x, a, b):
return a * (1 - np.exp(-b * x))
def monod(x, a, b, c, alp):
return c + (a * np.exp(b) * x) / (alp + x)
tic = time.perf_counter()
Vconvfact = -EL
vth = vtm/Vconvfact
vrm = vres/Vconvfact
t0_val = 0
vini_neg = EL
ts = np.inf
dt = d_dt/time_scale
print('dt= ', dt)
init_sign = 0
ref_t = 2
t0_val = 0
psi1 = ((-4)*bet+((1+delta1)**2))**(0.5)
Idep_ini = 0
Iadap_ini = 0
out = []
t_out = []
t_final = t0_val+dt
v_ini = -1
vini_prec = v_ini
mul = 15
v_star_min = -v_min/EL
alpha_neg = minCurr/sc
f = open(tSpikeOutputFileName, 'w')
i = 0
soglia_sign = 10
Ide = []
Iada = []
Ide2 = []
Iada2 = []
tetalist = []
t_spk = -3*d_dt
afirst = 0
meancorlastis = 0
stdcorlastis = 0
sis = 0
firstSpikeFlag = False
counter = 0
while(t_final*time_scale < sim_lenght):
if (t_final-init_sign)*time_scale >= tagliorette(cor[i],retteParParsed):
if cor[i] > ith:
if cor[i-1] < ith or i == 0:
init_sign = t_final
Idep_ini = max(Idep_ini_vr,cost_idep_ini*(cor[i]-ith))
Iadap_ini = 0
v_ini = ((EL + (1 - np.exp(-(zeta*1000*cor[i] - rho*1000*ith)/1000) )*(vtm - EL))/(-EL))
v_ini = V(t_final, delta1, psi1,
cor[i]/sc, bet, Iadap_ini, Idep_ini, t0_val, v_ini)
Iadap_ini = Iadap(
t_final, delta1, psi1, cor[i] / sc, bet, Iadap_ini, Idep_ini, t0_val, v_ini)
Idep_ini = Idep(t_final, bet, Idep_ini, t0_val)
if cor[i-1] == 0:
v_ini = vini_prec
else:
if cor[i]<ith and cor[i]>0:
v_ini = ((EL + (1 - np.exp(-(eta*1000*cor[i])/1000) )*(vtm - EL))/(-EL))
elif cor[i]<=0:
v_ini = ((EL + (1 - np.exp(-(csi*1000*cor[i])/1000) )*(vtm - EL))/(-EL))
else:
v_ini = ((EL + (1 - np.exp(-(zeta*1000*cor[i] - rho*1000*ith)/1000) )*(vtm - EL))/(-EL))
vini_prec = v_ini
out.append(v_ini)
t_out.append(t_final)
Iada.append(Iadap_ini)
Ide.append(Idep_ini)
else:
vini_prec = v_ini
if (cor[i] < ith and cor[i] >= 0) or i == 0:
v_ini = ((EL + (1 - np.exp(-(eta*1000*cor[i])/1000) )*(vtm - EL))/(-EL))
else:
if cor[i] < cor[i-1] and cor[i] > 0 and (t_spk+2*d_dt) < t_final*time_scale:
print('teta')
teta = (out[i-1]/(cor[i-1] / sc))*(1/dt-delta1) - \
(out[i-2]/((cor[i-1] / sc)*dt))-delta1/(cor[i-1] / sc)-1
if teta < 0:
teta = 0
Idep_ini = Iadap_ini + teta * (cor[i] / sc) / bet
tetalist.append(teta)
v_ini = V(t=t_final, delta=delta1, Psi=psi1,
alpha=cor[i]/sc, beta=bet, IaA0=Iadap_ini, IdA0=Idep_ini, t0=t0_val, V0=v_ini)
Iadap_ini = Iadap(t=t_final, delta=delta1, Psi=psi1,
alpha=cor[i]/sc, beta=bet, IaA0=Iadap_ini, IdA0=Idep_ini, t0=t0_val, V0=v_ini)
Idep_ini = Idep(t=t_final, beta=bet,
IdA0=Idep_ini, t0=t0_val)
else:
if cor[i] > 0:
v_ini = V(t=t_final, delta=delta1, Psi=psi1,
alpha=cor[i]/sc, beta=bet, IaA0=Iadap_ini, IdA0=Idep_ini, t0=t0_val, V0=v_ini)
Iadap_ini = Iadap(t=t_final, delta=delta1, Psi=psi1,
alpha=cor[i]/sc, beta=bet, IaA0=Iadap_ini, IdA0=Idep_ini, t0=t0_val, V0=v_ini)
Idep_ini = Idep(t=t_final, beta=bet,
IdA0=Idep_ini, t0=t0_val)
if cor[i-1] != cor[i] and (cor[i] < 0 and cor[i] > minCurr):
v_ini = vini_prec
if cor[i] < 0 and cor[i] > minCurr:
v_ini = ((EL + (1 - np.exp(-(csi*1000*cor[i])/1000) )*(vtm - EL))/(-EL))
if cor[i-1] != cor[i] and cor[i] <= minCurr:
Iadap_ini = -v_min/EL + 1
Idep_ini = 0
v_ini = ((EL + (1 - np.exp(-(csi*1000*cor[i])/1000) )*(vtm - EL))/(-EL))
if cor[i] <= minCurr:
v_ini=v_star_min
if v_ini*Vconvfact < v_min:
v_ini = v_min/Vconvfact
Iadap_ini = 0
out.append(v_ini)
t_out.append(t_final)
Iada.append(Iadap_ini)
Ide.append(Idep_ini)
if cor[i] > ith:
if cor[i-1] < ith:
init_sign = t_final
Idep_ini = max(Idep_ini_vr,cost_idep_ini*(cor[i]-ith))
Iadap_ini = 0
v_ini = ((EL + (1 - np.exp(-(zeta*1000*cor[i] - rho*1000*ith)/1000) )*(vtm - EL))/(-EL))
if cor[i-1]<1e-11:
v_ini = -1
if v_ini > vth:
t_spk = t_final*time_scale
f.write(str(round(t_spk, 3)) + ' \n')
v_ini = vrm
print('***spike***')
print('t ', t_final*time_scale, 'val_ist V', v_ini * Vconvfact, 'adap',
Iadap_ini, 'idep', Idep_ini, 't_ini', init_sign)
print('************')
if cor[i] < istim_min_spikinig_exp or cor[i] > istim_max_spikinig_exp:
c_aux = c
Iadap_ini = monod((t_final-init_sign) *
time_scale, a, b*cor[i]/1000, c_aux, alp)
else:
Iadap_ini = monod((t_final-init_sign) * time_scale, a, b*cor[i]/1000, c, alp)
if Iadap_ini<0:
print('monod negativa')
paramL = Iadap_ini
if a > 0:
c_aux = c - paramL
else:
c_aux = -a*np.exp(b*cor[i]/1000)
Iadap_ini = monod((t_final-init_sign) * time_scale, a, b*cor[i]/1000, c_aux, alp)
if cor[i] < ith:
Idep_ini = 0
Iadap_ini = 0
else:
Idep_ini = Idep_ini_vr
for k in range(int(ref_t / d_dt)):
out.append(v_ini)
t_final = t_final + dt
t_out.append(t_final)
Iada.append(Iadap_ini)
Ide.append(Idep_ini)
i = i + 1
vini_prec = v_ini
i = i + 1
t0_val = t_final
t_final = t0_val+dt
file = open(voltageOutputFileName, "w")
for i in range(len(t_out)):
file.write(str(t_out[i]*time_scale) + " " + str(out[i]*Vconvfact) + "\n")
file.close()
print(t_final)
toc = time.perf_counter()
print(f"time: {toc - tic:0.4f} seconds")
plt.show()
f.close()
'''
MAIN
'''
# user input -------------------------------------------------------------------
currentFileName = 'trace_l23-06-13.res.6-tt6clu6_0.65+1_prova20122022.txt'
tSpikeOutputFileName = '6_065_t_spk_simulated.txt'
voltageOutputFileName = '6_065_voltage_simulated.txt'
# Neuron parameters
EL = -72.5
vres = -65
vtm = -52
Cm = 2047.4164432004916
ith = 300.9987901902274
tao_m = 3310.462136574088
sc = 4526.328798037026
bet = 0.24522251335200956
delta1 = 0.009906254244852036
cost_idep_ini = 0.017625482908326662
Idep_ini_vr = 1.0122905259090516
psi1 = 0.1975362978159442
a=14.2787
b=-2.10966
c=0.0608809
alp=225.491
istim_min_spikinig_exp=400
istim_max_spikinig_exp=1000
sim_lenght = 600000.1#milliseconds
rettePar=[np.inf,-0.028,76.86,700,0.68,-190]
# equilibrium parameters
v_min = -90
minCurr = -185
# Neuron
zeta = 3.5e-3
eta = 2.5e-3
rho = 1e-3
csi = 3.5e-3
# -------------------------------------------------------------------
neuronParameters=[EL,vres,vtm,Cm,ith,tao_m,sc,bet,delta1,cost_idep_ini,Idep_ini_vr,psi1,a,b,c,alp,istim_min_spikinig_exp,istim_max_spikinig_exp,sim_lenght,rettePar]
equilibriumParameters = [v_min,minCurr,zeta,eta,rho,csi]
AGLIFsynaptic(neuronParameters,equilibriumParameters,currentFileName,tSpikeOutputFileName,voltageOutputFileName)