# ----------------------------------------------------------
# Main simulation control
#
# Naoki Hiratani, N.Hiratani@gmail.com
#
# Based on the model by Dr. Tiago Branco (Smith et al., 2013)
#
# ----------------------------------------------------------
import numpy as np
import pickle
import time
import libcell as lb
#import saveClass as sc
from math import *
import matplotlib.pyplot as plt
import random as pyrnd
import sys
#----------------------------------------------------------------------------
# Functions and Classes
#calculation of maximum membrane potential difference
def calc_vks(v,Nin):
vkmin = min(v)
vkmaxs = []
vks = []
for tidx in range(2,len(v)):
if (v[tidx-2] - v[tidx-1])*(v[tidx-1] - v[tidx]) < 0.0:
if v[tidx-1] > vkmin + 0.01:
vks.append( v[tidx-1] - v[-1] )
if len(vks) != Nin:
print "the program failed to estimate vks"
return vks
#calculation of unit EPSPs
def estimate_vkappas(model):
locs = lb.calc_locs_kappa(model)
lb.add_Estims(model,locs)
lb.add_AMPAsyns(model,locs,model.gmax)
model.tstop = 100
strvkp = "data/spf4d4_vkappas_RA" + str(model.RA) + "_gm" + str(model.gmax) + "_dm_" + str(model.dendmodel) + ".txt"
fvkp = open(strvkp,'w')
vkappas = []
for lidx in range(len(locs)):
for l2idx in range(len(locs)):
model.Estimlist[l2idx].number = 0
model.Estimlist[lidx].number = 1
ttmp, vtmp = lb.simulate(model)
dvtmp = max(vtmp) - vtmp[-1]
print lidx,dvtmp
vkappas.append( dvtmp )
fvkp.write( str(locs[lidx][0]) + " " + str(locs[lidx][1]) + " " + str(dvtmp) + "\n" )
fvkp.flush()
plt.hist(vkappas)
plt.show()
def estimate_vks(model, vks):
model.sinterval = 1000.0
for i in range(model.Nin):
model.Estimlist[i].interval = model.sinterval
model.Estimlist[i].number = 1
model.Estimlist[i].start = 500 + i*model.sinterval
model.Estimlist[i].noise = 0
model.tstop = model.Nin*model.sinterval + 1000
ttmp, vtmp = lb.simulate(model)
return calc_vks(vtmp,model.Nin)
def pickup_vks(model,locs): #pick up the value of vk from the file 'strvkp'
strvkp = "data/nrn_simul_vkappas_RA" + str(model.RA) + "_gm" + str(model.gmax) + "_dm_" + str(model.dendmodel) + ".txt"
model.dendidxs = []
model.dendpos = []
model.vkps = []
for line in open(strvkp,'r'):
vktmp = line.split(" ")
model.dendidxs.append( int(vktmp[0]) )
model.dendpos.append( float(vktmp[1]) )
model.vkps.append( float(vktmp[2].strip("\n")) )
vks = []
for nidx in range(len(locs)):
lidx = locs[nidx][0]
ptmp = locs[nidx][1]
dptmp = 1.0
vks.append(0.0)
for n2idx in range(len(model.dendidxs)):
if model.dendidxs[n2idx] == lidx:
if abs( ptmp - model.dendpos[n2idx] ) <= dptmp:
vks[nidx] = model.vkps[n2idx]
dptmp = abs( ptmp - model.dendpos[n2idx] )
if vks[nidx] == 0.0:
print "an error occurred at nidx = " + str(nidx) + ", while reading the vkappa file."
#print vks
return vks
def pick_a_vk(model,loc):
vktmp = 0.0
lidx = loc[0]
ptmp = loc[1]
dptmp = 1.0
for nidx in range(len(model.dendidxs)):
if model.dendidxs[nidx] == lidx:
if abs( ptmp - model.dendpos[nidx] ) <= dptmp:
vktmp = model.vkps[nidx]
dptmp = abs( ptmp - model.dendpos[nidx] )
return vktmp
def estimate_dks(model,locs): #distance from the soma to the synapses
dks = []
for lidx in range(len(locs)):
loc = locs[lidx]
dks.append( model.dists[loc[0]] + loc[1]*model.dends[loc[0]].L )
return dks
def calc_pvks(vks): #distribution of vk
pvks = []
vkmin = min(vks)
vkmax = max(vks)
dvk = (vkmax - vkmin)/10.0
for i in range(len(vks)):
pvks.append(0.0)
for i2 in range(len(vks)):
if vks[i]-0.5*dvk < vks[i2] and vks[i2] < vks[i]+0.5*dvk:
pvks[i] += 1.0/float(len(vks)*dvk)
return pvks
def initialize_ukgks(model,pvks):
uks = []
gks = []
Zuks = []
for j in range(model.Min):
Zuks.append(0.0)
for i in range(model.Nin):
uks.append(1.0/float(pvks[i]))
Zuks[ model.preidx[i] ] += uks[i]
gks.append(1.0)
for i in range(model.Nin):
uks[i] = uks[i]/Zuks[ model.preidx[i] ]
return uks,gks
#update of spine sizes with near-optimal learning
def update_uks(model, vks, uks, pre_spikes, rewiring, dropping):
gamma = model.gamma
uk_min = model.uk_min; uk_max = model.uk_max
vkmin = model.vkmin
min_rate = model.min_rate
prev_uks = uks
connec_exists = []
for j in range(model.Min):
connec_exists.append(0)
for i in range(model.Nin):
jidx = model.preidx[i];
if uks[i] > 0.0:
connec_exists[jidx] += 1
Zuks = np.zeros((model.Min))
for i in range(model.Nin):
jidx = model.preidx[i]
if connec_exists[jidx] > 0:
tmp_rate = min_rate*exp( gamma*vks[i] )
uks[i] *= (tmp_rate**pre_spikes[i])*exp(-tmp_rate)/factorial( pre_spikes[i] )
Zuks[jidx] += uks[i]
for i in range(model.Nin):
jidx = model.preidx[i]
if connec_exists[jidx] > 0:
if Zuks[jidx] > 0.0:
uks[i] = uks[i]/Zuks[jidx]
else:
uks[i] = prev_uks[i]
if uks[i] > 0.0:
print model.preidx[i], i, vks[i], uks[i]
if prev_uks[i] > 0.0 and uks[i] == 0.0:
uks[i] = 0.1*uk_min
if (not rewiring) and (not dropping) and uks[i] < uk_min:
uks[i] = uk_min
if uks[i] > uk_max:
uks[i] = uk_max
return uks
#update of spine sizes without pre-neuron specific normalization
def update_uks_apr(model, vks, uks, prespikes, rewiring, dropping):
Nin = model.Nin; Min = model.Min; Kin = model.Kin
gamma = model.gamma
uk_min = model.uk_min; uk_max = model.uk_max
vkmin = model.vkmin; vo = model.vo
min_rate = model.min_rate
duks = []
Zdtot = 0.0
for i in range(Nin):
jidx = model.preidx[i]
rhojk = min_rate*exp(gamma*vks[i])
vojk = (1.0 - uks[i])*vo + uks[i]*vks[i]
rhoo = min_rate*exp(gamma*vojk)
duks.append( ( (rhojk/rhoo)**prespikes[i] )*exp( -(rhojk-rhoo) ) )
Zdtot += log(duks[i])/float(Nin)
epZdtot = exp(Zdtot)
for i in range(Nin):
uks[i] = uks[i]*duks[i]/epZdtot
if uks[i] > 0.5*uk_max:
uks[i] = 0.5*uk_max
if (not rewiring) and (not dropping) and uks[i] < uk_min:
uks[i] = uk_min
return uks
#calculation of EPSP area
def calc_integ_v(vtmp,drt):
integ_v = 0.0
vmin = vtmp[400] #min(vtmp) # membrane potential at t=10ms
for tidx in range(len(vtmp)):
integ_v += (vtmp[tidx] - vmin)*drt
return integ_v
def evaluate_perf(model,prespikes,ttmp,vtmp):
epsilon = model.epsilon
opt_pv = 0.0; opt_nv = 0.0;
for j in range(model.Min):
opt_pv += prespikes[j]*(model.log_p_rates[j] + log(epsilon))
opt_nv += prespikes[j]*(model.log_n_rates[j] + log(epsilon))
drt = ttmp[1] - ttmp[0]
integ_v = calc_integ_v(vtmp,drt)
depol_v = max(vtmp) - vtmp[400]
total_prespikes = 0
for prespike in prespikes:
total_prespikes += prespike
perftmp = [opt_pv, opt_nv, integ_v, depol_v, total_prespikes]
return perftmp
def is_match(t,trs):
tof = False
for tr in trs:
if t == tr:
tof = True
return tof
def simul(Kin, gmax, gI, uk_min, sd_bias, release_prob, ik):
model = lb.L23() # Layer 2/3 model is selected
dendActive = True # with passive/active dendrrite
opt_update = True #whether update rule is optimal or approximate
rewiring = True #if the model includes rewiring or not
dropping = True #if the model includes uncompensated elimination
membrane_recording = True #if membrane dynamics is recorded or not (generate a heavy file if true)
prespike_recording = False #if the prespikes are recorded or not
# active channels are distributed over both axon and dendrite
lb.init_active(model, axon=True, soma=True, dend=dendActive, dendNa=True, dendCa=True)
lb.init_params(model, Kin, gmax, gI, uk_min, sd_bias, release_prob)
model.dendmodel = 'active' if dendActive else 'passive'
model.learningrule = 'opt' if opt_update else 'apr'
if rewiring or dropping:
model.connections = ''
if rewiring:
model.connections = model.connections + 'rewired_'
if dropping:
model.connections = model.connections + 'dropped_'
else:
model.connections = 'fixed'
#CAUTION: 'estimate_vkappas' resets vkappa
#estimate_vkappas(model)
#locs = lb.calc_locs(model) #uniformly locate synapses
locs = lb.calc_locs_random(model) # randomly locate synapses
lb.allocate_syns(model)
#lb.plot_morphology(model,locs)
lb.add_Estims(model,locs)
lb.add_AMPAsyns(model, locs, model.gmax)
inh_locs = lb.calc_inh_locs_uniform(model)
lb.add_Istims(model, inh_locs)
lb.add_GABAsyns(model, inh_locs, model.gI)
dks = estimate_dks(model,locs) #estimate the distances from the soma
vks = pickup_vks(model,locs) #readout values of vks
mean_pre_rates = np.full((model.Nin), model.spn_rate)
pvks = calc_pvks(vks) #prior distribution of vks
uks, gks = initialize_ukgks(model,pvks)
lb.generate_pre_tuning(model)
lb.generate_input_rates(model)
model.vkmax = max(vks); model.vkmin = min(vks)
model.gamma = max(model.log_p_rates)/model.vkmax
model.vo = 1.5*model.vkmin#
fstrtmp = 'K' + str(model.Kin) + "_RA" + str(model.RA) + '_gm' + str(model.gmax) + '_gi' + str(model.gI) + \
'_dm_' + str(model.dendmodel) + '_lr_' + str(model.learningrule) + '_cw_' + str(model.connections) \
+ '_ukm_' + str(uk_min) + '_sdb_' + str(sd_bias) + '_pr_' + str(release_prob) + '_ik' + str(ik) + '.txt'
strperf = 'data/nrn_simul_perf_' + fstrtmp
fperf = open(strperf,'w')
strduvk = 'data/nrn_simul_duvk_' + fstrtmp
fduvk = open(strduvk,'w')
strtune = 'data/nrn_simul_tune_' + fstrtmp
ftune = open(strtune,'w')
for j in range(model.Min):
ftune.write(str(j) + ' ' + str(model.phis[j]) + ' ' + str(model.thetas[j]) + ' ' + str(model.rfdists[j]) + ' ' + str(model.log_p_rates[j]/model.gamma) + ' ' + str(model.log_n_rates[j]/model.gamma) + '\n')
ftune.flush()
if prespike_recording:
strprespike = 'data/nrn_simul_prespike_' + fstrtmp
fprespike = open(strprespike,'w')
model.tstop = 100
trs = [0,100]
perfs = []
for q in range(3):
perfs.append([])
for t in range(model.trials):
if is_match(t,trs):#evaluation
if membrane_recording:
if t == trs[0]:
strmbps = 'data/nrn_simul_mbps_' + fstrtmp
fmbps = open(strmbps,'w')
elif t == trs[-1]:
strmbpf = 'data/nrn_simul_mbpf_' + fstrtmp
fmbpf = open(strmbpf,'w')
for i in range(model.Nin):
fduvk.write(str(locs[i][0]) + " " + str(locs[i][1]) + " " + str(model.preidx[i]) + " " + str(dks[i]) + " " + str(uks[i]) + " " + str(vks[i]) + "\n")
dist_bt_syns = lb.calc_dist_bt_syns(model, locs)
for didx in range(len(dist_bt_syns)):
fduvk.write( str(dist_bt_syns[didx]) + ' ' )
fduvk.write( '\n' ); fduvk.flush()
#test phase
for t2 in range(200):
true_pre_spikes = []; pre_spikes = []
if t2 < 100:
true_pre_spikes, pre_spikes = lb.generate_prespikes(model, model.p_rates)
else:
true_pre_spikes, pre_spikes = lb.generate_prespikes(model, model.n_rates)
if prespike_recording:
for j in range(model.Min):
fprespike.write( str(t) + " " + str(t2) + " " + str(j) + " " + str(true_pre_spikes[j]) + "\n" )
fprespike.flush()
lb.generate_inputs(model, uks, pre_spikes)
lb.generate_inh_inputs(model, true_pre_spikes)
ttmp, vtmp = lb.simulate(model)
perftmp = evaluate_perf(model,true_pre_spikes,ttmp,vtmp)
strtmp = str(t) + ' ' + str(perftmp[0]) + ' ' + str(perftmp[1]) + ' ' + str(perftmp[2]) + ' ' + str(perftmp[3]) + ' ' + str(perftmp[4]) + '\n'
fperf.write(strtmp)
if membrane_recording:
if t == trs[0]:
for t3 in range(len(ttmp)):
strtmp = str(t2) + ' ' + str(ttmp[t3]) + ' ' + str(vtmp[t3]) + '\n'
fmbps.write( strtmp )
fmbps.flush()
if t == trs[-1]:
for t3 in range(len(ttmp)):
strtmp = str(t2) + ' ' + str(ttmp[t3]) + ' ' + str(vtmp[t3]) + '\n'
fmbpf.write( strtmp )
fmbpf.flush()
#training phase
true_pre_spikes, pre_spikes = lb.generate_prespikes(model, model.p_rates)
mean_pre_rates = np.multiply(1.0 - 1.0/model.taumr, mean_pre_rates) + np.multiply(1.0/model.taumr, pre_spikes)
if opt_update:
uks = update_uks(model, vks, uks, pre_spikes, rewiring, dropping)
else:
uks = update_uks_apr(model, vks, uks, pre_spikes, rewiring, dropping)
if prespike_recording:
for j in range(model.Min):
fprespike.write( str(t) + " " + str(j) + " " + str(true_pre_spikes[j]) + "\n" )
fprespike.flush()
if rewiring: #rewiring of connections with small spine sizes
rwcnt = 0
loctmp = []
for i in range(model.Nin):
if uks[i] < model.ukthreshold and uks[i] != 0.0:
if np.random.random() < model.rewiring_freq:
loctmp = lb.restricted_resampling(model, locs, i)
locs[i][0] = loctmp[0]; locs[i][1] = loctmp[1]
vks[i] = pick_a_vk(model,locs[i])
uks[i] = model.ukinit
dks[i] = model.dists[locs[i][0]] + locs[i][1]*model.dends[locs[i][0]].L
lb.rewire_synapse(model,locs[i],i)
rwcnt += 1
if dropping: #elimination of inactive connections
for i in range(model.Nin):
if mean_pre_rates[i] < model.pre_rates_th and np.random.random() < model.rewiring_freq:
uks[i] = 0.0
for i in range(model.Nin):
if uks[i] != 0.0 and uks[i] < model.ukthreshold:
uks[i] = model.uk_min
if __name__ == "__main__":
#def main(args=None):
#"""Main"""
param = sys.argv
Kin = int(param[1]) #total redundancy of synapses (Kin=5)
gmax = float(param[2]) #standard conductance [muS] (gmax=0.0025)
gI = float(param[3]) #inhibitory conductance [muS] (gI = 0.00075)
uk_min = float(param[4]) #threshold of rewiring (ukth=0.001)
sd_bias = float(param[5]) #bias in synaptic distribution (sd_bias=1.0)
release_prob = float(param[6]) #Release probability of synaptic vesicle (release_prob = 1.0)
ik = 0
simul(Kin, gmax, gI, uk_min, sd_bias, release_prob, ik)