from neuron import h,rxd
from pylab import *
import mytools
import scipy.io
import re
stim_T = 300.0
Nstims = 2
decay = 80
axonPasLen = 80
axonActLen = 20
Vesmax = 0.1 #Concentration of vesicles (mM)
Fmax = 0.001 #Concentration of fusion factor (mM)
initfile = 'None'
tolstochange = 'Gi,GiaGDP,Gibg'.split(',')
tolschange = [float(x) for x in '2,2,2'.split(',')]
gaba_input_dur = 3.0
blockedStr = 'None'
blockeds = []
block_factors = []
isyntype = 1 #1 for E->E connection (OnlyExp1_RGS used), 2 for I->E connection (OnlyExp2_RGS used)
spineLen = 0.8
spineDiam = 0.4
doInterpolate = 1
doDraw = 0
fileOut = 'None'
stimAmp = 0.3
kcoeff_f = 1.0
kcoeff_b = 1.0
kcoeff_2f = 1.0
kcoeff_2b = 1.0
kcoeff_3 = 1.0
kcoeff_4 = 1.0
gnabar = 3.9 #L5PCbiophysics4.hoc
gkbar = 0.077 #L5PCbiophysics4.hoc
gl = 0.00003 #L5PCbiophysics4.hoc
el = -90 #L5PCbiophysics4.hoc
gCaNtypebar_single = 2.7e-6 #uS (https://www.nature.com/articles/nn.2657) #13e-6 #uS (https://www.ncbi.nlm.nih.gov/books/NBK6181/)
minCa = 50e-6
params_all = ['107.034,514.651,167.633,11.026,66.931,3.009,1.904,2.241,1507.740','91.805,846.811,67.856,9.858,71.217,5.000,1.817,2.170,3000.000','133.092,528.522,185.581,8.574,25.868,3.186,2.125,2.729,2815.426','108.557,729.183,292.710,2.473,16.923,4.406,2.127,2.973,696.729','111.531,708.687,165.616,10.797,93.218,2.903,1.728,2.235,2797.567','109.869,469.365,114.137,59.053,73.263,3.573,2.385,3.028,2852.993','84.387,57.393,126.013,4.253,57.544,5.000,1.722,2.122,929.204','126.993,135.047,348.945,3.733,18.257,3.800,2.647,3.237,2804.842','105.438,642.913,142.107,19.602,40.482,3.599,2.391,2.974,1386.488','91.403,459.570,97.658,3.443,65.833,5.000,1.540,1.942,2447.114','71.410,380.140,143.521,2.991,45.796,10.000,2.044,2.836,2126.204','146.055,702.846,310.725,3.041,13.076,3.334,2.351,3.166,2112.405','90.417,993.718,98.403,6.879,95.287,3.944,1.572,2.204,1413.654','119.681,591.014,91.960,89.687,88.948,3.453,2.300,2.755,2950.763','123.170,616.560,117.913,74.842,46.625,3.952,2.609,3.283,3000.000']
params_str = 'p0'
if len(sys.argv) > 1:
gnabar = float(sys.argv[1])
if len(sys.argv) > 2:
gkbar = float(sys.argv[2])
if len(sys.argv) > 3:
gl = float(sys.argv[3])
if len(sys.argv) > 4:
el = float(sys.argv[4])
if len(sys.argv) > 5:
decay = float(sys.argv[5])
if len(sys.argv) > 6:
minCa = float(sys.argv[6])
if len(sys.argv) > 7:
params_str = sys.argv[7]
if len(sys.argv) > 8:
blockedStr = sys.argv[8]
if len(sys.argv) > 9:
Nstims = int(float(sys.argv[9]))
if len(sys.argv) > 10:
stim_T = float(sys.argv[10])
if len(sys.argv) > 11:
axonPasLen = int(float(sys.argv[11]))
if len(sys.argv) > 12:
axonActLen = int(float(sys.argv[12]))
if len(sys.argv) > 13:
isyntype = int(float(sys.argv[13]))
if len(sys.argv) > 14:
stimAmp = float(sys.argv[14])
if len(sys.argv) > 15:
kcoeff_f = float(sys.argv[15])
if len(sys.argv) > 16:
kcoeff_b = float(sys.argv[16])
if len(sys.argv) > 17:
kcoeff_2f = float(sys.argv[17])
if len(sys.argv) > 18:
kcoeff_2b = float(sys.argv[18])
if len(sys.argv) > 19:
kcoeff_3 = float(sys.argv[19])
if len(sys.argv) > 20:
kcoeff_4 = float(sys.argv[20])
if len(sys.argv) > 21:
Fmax = float(sys.argv[21])
if len(sys.argv) > 22:
Vesmax = float(sys.argv[22])
if len(sys.argv) > 23:
doInterpolate = int(float(sys.argv[23]))
if len(sys.argv) > 24:
doDraw = int(float(sys.argv[24]))
if len(sys.argv) > 25:
fileOut = sys.argv[25]
if params_str[0] == 'p':
iparam = int(params_str[1:])
params = {'k[0]': float(params_all[iparam].split(',')[0]), 'k[1,9]': float(params_all[iparam].split(',')[1]), 'k[3,4,5,6,7,8]': float(params_all[iparam].split(',')[2]), 'k[15,17,19,21,23]': float(params_all[iparam].split(',')[3]), 'k[16,18,20,22,24]': float(params_all[iparam].split(',')[4]), 'OnlyExp0_RGS': float(params_all[iparam].split(',')[5]), 'OnlyExp1_RGS': float(params_all[iparam].split(',')[6]), 'OnlyExp2_RGS': float(params_all[iparam].split(',')[7]), 'gaba_flux': float(params_all[iparam].split(',')[8])}
else:
params = {'k[0]': float(params_str.split(',')[0]), 'k[1,9]': float(params_str.split(',')[1]), 'k[3,4,5,6,7,8]': float(params_str.split(',')[2]), 'k[15,17,19,21,23]': float(params_str.split(',')[3]), 'k[16,18,20,22,24]': float(params_str.split(',')[4]), 'OnlyExp0_RGS': float(params_str.split(',')[5]), 'OnlyExp1_RGS': float(params_str.split(',')[6]), 'OnlyExp2_RGS': float(params_str.split(',')[7]), 'gaba_flux': float(params_str.split(',')[8])}
gaba_input_flux = params['gaba_flux']
styles = ['-','--']
widths = [0.3,0.6]
tstop = 3000
if 1000+stim_T*Nstims > 3000:
tstop = 1000+stim_T*Nstims
#Create the axon and the presynaptic terminal
h("""
load_file("stdlib.hoc")
load_file("stdrun.hoc")
objref cvode
cvode = new CVode()
cvode.active(1)
cvode.atol(0.0000005)
tstop = """+str(tstop)+"""
create axon[4]
create spine
connect axon[1](0), axon[0](1)
connect axon[2](0), axon[1](1)
connect axon[3](0), axon[2](1)
connect spine(0), axon[1](1)
for (i=0;i<4;i+=3) {
axon[i] L = """+str(axonPasLen)+"""
axon[i] diam = 1.0
nseg = 5
}
for (i=1;i<3;i+=1) {
axon[i] L = """+str(axonActLen)+"""
axon[i] diam = 1.0
axon[i] insert NaTa_t
axon[i] insert K_Tst
axon[i] insert pas
axon[i] insert Ih
axon[i] insert Im
axon[i] insert CaDynamics_E2
axon[i] insert Ca_LVAst
axon[i] insert Ca_HVA
axon[i] insert SKv3_1
axon[i] insert SK_E2
axon[i] insert K_Pst
axon[i] insert Nap_Et2
axon[i] gNaTa_tbar_NaTa_t = """+str(gnabar)+"""
axon[i] gK_Tstbar_K_Tst = """+str(gkbar)+"""
axon[i] g_pas = """+str(gl)+"""
axon[i] e_pas = """+str(el)+"""
axon[i] gIhbar_Ih = 0.0001
axon[i] gImbar_Im = 0.013322
axon[i] decay_CaDynamics_E2 = 277.300774
axon[i] gamma_CaDynamics_E2 = 0.000525
axon[i] gCa_LVAstbar_Ca_LVAst = 0.000813
axon[i] gCa_HVAbar_Ca_HVA = 0.000222
axon[i] gSKv3_1bar_SKv3_1 = 0.473799
axon[i] gSK_E2bar_SK_E2 = 0.000047
axon[i] gK_Pstbar_K_Pst = 0.188851
axon[i] gNap_Et2bar_Nap_Et2 = 0.005834
nseg = 5
}
spine L = """+str(spineLen)+"""
spine.diam = """+str(spineDiam)+"""
spine insert CaNtypeWillingReluctant
spine nseg = 1
""")
cyt = rxd.Region([h.spine], name='cyt', nrn_region='i')
my_volume = spineLen * pi * (spineDiam/2)**2 * 1e-15 #in litres; multiplied by 1e-18 to convert um3->m3 and by 1e3 for m3->l
ks = [1.0]*31
ks[0] = 0.0005 # gaba --> gabaOut (forward)
ks[1] = 5.555000000000001 # gaba + GABABR <-> gabaGABABR (forward)
ks[2] = 0.005 # gaba + GABABR <-> gabaGABABR (backward)
ks[3] = 150.0 # gabaGABABR + Gi <-> gabaGABABRGi (forward)
ks[4] = 0.00025 # gabaGABABR + Gi <-> gabaGABABRGi (backward)
ks[5] = 0.000125 # gabaGABABRGi --> gabaGABABRGibg + GiaGTP (forward)
ks[6] = 0.001 # gabaGABABRGibg --> gabaGABABR + Gibg (forward)
ks[7] = 75.0 # GABABR + Gi <-> GABABRGi (forward)
ks[8] = 0.000125 # GABABR + Gi <-> GABABRGi (backward)
ks[9] = 5.555000000000001 # gaba + GABABRGi <-> gabaGABABRGi (forward)
ks[10] = 0.005 # gaba + GABABRGi <-> gabaGABABRGi (backward)
ks[11] = 2.0 # GiaGTP + RGS <-> GiaGTPRGS (forward)
ks[12] = 0.002 # GiaGTP + RGS <-> GiaGTPRGS (backward)
ks[13] = 0.03 # GiaGTPRGS --> GiaGDP + RGS (forward)
ks[14] = 1250.0 # GiaGDP + Gibg --> Gi (forward)
ks[15] = 14.0 # GIRK + Gibg <-> GIRKGibg (forward)
ks[16] = 0.001 # GIRK + Gibg <-> GIRKGibg (backward)
ks[17] = 14.0 # GIRKGibg + Gibg <-> GIRKGibg2 (forward)
ks[18] = 0.001 # GIRKGibg + Gibg <-> GIRKGibg2 (backward)
ks[19] = 14.0 # GIRKGibg2 + Gibg <-> GIRKGibg3 (forward)
ks[20] = 0.001 # GIRKGibg2 + Gibg <-> GIRKGibg3 (backward)
ks[21] = 14.0 # GIRKGibg3 + Gibg <-> GIRKGibg4 (forward)
ks[22] = 0.001 # GIRKGibg3 + Gibg <-> GIRKGibg4 (backward)
ks[23] = 14.0 # VGCC + Gibg <-> VGCCGibg (forward)
ks[24] = 0.001 # VGCC + Gibg <-> VGCCGibg (backward)
ks[25] = 1e16*kcoeff_f # 4*ca + Ffactor <-> Ffactorca4 (forward)
ks[26] = 0.1*kcoeff_b # 4*ca + Ffactor <-> Ffactorca4 (backward)
ks[27] = 1000*kcoeff_2f # Ffactorca4 + Ves <-> Ffactorca4Ves (forward)
ks[28] = 0.1*kcoeff_2b # Ffactorca4 + Ves <-> Ffactorca4Ves (backward)
ks[29] = 4*kcoeff_3 # Ffactorca4Ves --> 10000*NT + Ffactor + Ves + 4*ca
ks[30] = 10*kcoeff_4 # NT -->
initvalues = [0.0, 0.0, 0.00039999999999999996, 0.0, 0.0026, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.001, 9.999999999999999e-05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, Fmax, 0.0, Vesmax, 0.0, 0.0]
#Use the fitted values
initvalues[10] = initvalues[10]*params['OnlyExp'+str(isyntype)+'_RGS'] #Presynaptic, E->E or I->E connection depending on type
ks[0] = ks[0]*params['k[0]']
for i in [1,9]:
ks[i] = ks[i]*params['k[1,9]']
for i in [2,10]:
ks[i] = ks[i]*5.555*params['k[1,9]']*0.00011/0.005
for i in [3,4,5,6,7,8]:
ks[i] = ks[i]*params['k[3,4,5,6,7,8]']
for i in [15,17,19,21,23]:
ks[i] = ks[i]*params['k[15,17,19,21,23]']
for i in [16,18,20,22,24]:
ks[i] = ks[i]*params['k[16,18,20,22,24]']
initvalues[13] = 0 # No GIRKs presynaptically
#Create the species and set their initial values. Set the tighter relative tolerances for Ca-dependent release
species = ['gaba', 'gabaOut', 'GABABR', 'gabaGABABR', 'Gi', 'GABABRGi', 'gabaGABABRGi', 'gabaGABABRGibg', 'GiaGTP', 'GiaGDP', 'RGS', 'GiaGTPRGS', 'Gibg', 'GIRK', 'VGCC', 'GIRKGibg', 'GIRKGibg2', 'GIRKGibg3', 'GIRKGibg4', 'VGCCGibg', 'ca', 'Ffactor', 'Ffactorca4', 'Ves', 'Ffactorca4Ves', 'NT']
print("my_volume = "+str(my_volume)+" l ?= "+str(h.spine.L*(h.spine.diam/2)**2*3.14159265358)+" um3")
if len(initfile) > 3 and initfile != 'None':
DATA_init = scipy.io.loadmat(initfile)
for ispec in range(0,len(species)):
initvalues[ispec] = DATA_init['DATA'][1+ispec,-1]
if species[ispec] not in DATA_init['headers'][1+ispec]:
print("Warning: mismatch in DATA_init, ispec = "+str(ispec))
tolscales = [1.0 for i in range(0,len(species))]
for ispec in [20,21,22,23,24]:
tolscales[ispec] = 0.001
if blockedStr.find('x') > -1:
blockeds = blockedStr[0:blockedStr.rfind('x')].split(',')
block_factors = [float(x) for x in blockedStr[blockedStr.rfind('x')+1:].split(',')]
print('blockeds = '+str(blockeds)+', block_factors = '+str(block_factors))
for ispec in range(0,len(species)):
for iblock in range(0,len(blockeds)):
if re.match(blockeds[iblock]+'$',species[ispec]):
initvalues[ispec] = block_factors[iblock]*initvalues[ispec]
print('Match found, ispec = '+str(ispec))
for itol in range(0,len(tolstochange)):
if tolstochange[itol] == species[ispec]:
tolscales[ispec] = tolscales[ispec]*10**(-tolschange[itol])
for iblock in range(0,len(blockeds)):
if blockeds[iblock] == 'gaba_flux':
gaba_input_flux = gaba_input_flux*block_factors[iblock]
gaba = rxd.Species(cyt, name='gaba', charge=0, initial=initvalues[0], atolscale=tolscales[0])
gabaOut = rxd.Species(cyt, name='gabaOut', charge=0, initial=initvalues[1], atolscale=tolscales[1])
GABABR = rxd.Species(cyt, name='GABABR', charge=0, initial=initvalues[2], atolscale=tolscales[2])
gabaGABABR = rxd.Species(cyt, name='gabaGABABR', charge=0, initial=initvalues[3], atolscale=tolscales[3])
Gi = rxd.Species(cyt, name='Gi', charge=0, initial=initvalues[4], atolscale=tolscales[4])
GABABRGi = rxd.Species(cyt, name='GABABRGi', charge=0, initial=initvalues[5], atolscale=tolscales[5])
gabaGABABRGi = rxd.Species(cyt, name='gabaGABABRGi', charge=0, initial=initvalues[6], atolscale=tolscales[6])
gabaGABABRGibg = rxd.Species(cyt, name='gabaGABABRGibg', charge=0, initial=initvalues[7], atolscale=tolscales[7])
GiaGTP = rxd.Species(cyt, name='GiaGTP', charge=0, initial=initvalues[8], atolscale=tolscales[8])
GiaGDP = rxd.Species(cyt, name='GiaGDP', charge=0, initial=initvalues[9], atolscale=tolscales[9])
RGS = rxd.Species(cyt, name='RGS', charge=0, initial=initvalues[10], atolscale=tolscales[10])
GiaGTPRGS = rxd.Species(cyt, name='GiaGTPRGS', charge=0, initial=initvalues[11], atolscale=tolscales[11])
Gibg = rxd.Species(cyt, name='Gibg', charge=0, initial=initvalues[12], atolscale=tolscales[12])
GIRK = rxd.Species(cyt, name='GIRK', charge=0, initial=initvalues[13], atolscale=tolscales[13])
VGCC = rxd.Species(cyt, name='VGCC', charge=0, initial=initvalues[14], atolscale=tolscales[14])
GIRKGibg = rxd.Species(cyt, name='GIRKGibg', charge=0, initial=initvalues[15], atolscale=tolscales[15])
GIRKGibg2 = rxd.Species(cyt, name='GIRKGibg2', charge=0, initial=initvalues[16], atolscale=tolscales[16])
GIRKGibg3 = rxd.Species(cyt, name='GIRKGibg3', charge=0, initial=initvalues[17], atolscale=tolscales[17])
GIRKGibg4 = rxd.Species(cyt, name='GIRKGibg4', charge=0, initial=initvalues[18], atolscale=tolscales[18])
VGCCGibg = rxd.Species(cyt, name='VGCCGibg', charge=0, initial=initvalues[19], atolscale=tolscales[19])
ca = rxd.Species(cyt, name='ca', charge=2, initial=initvalues[20], atolscale=tolscales[20])
Ffactor = rxd.Species(cyt, name='Ffactor', charge=0, initial=initvalues[21], atolscale=tolscales[21])
Ffactorca4 = rxd.Species(cyt, name='Ffactorca4', charge=0, initial=initvalues[22], atolscale=tolscales[22])
Ves = rxd.Species(cyt, name='Ves', charge=0, initial=initvalues[23], atolscale=tolscales[23])
Ffactorca4Ves = rxd.Species(cyt, name='Ffactorca4Ves', charge=0, initial=initvalues[24], atolscale=tolscales[24])
NT = rxd.Species(cyt, name='NT', charge=0, initial=initvalues[25], atolscale=tolscales[25])
gaba_flux_rate = rxd.Parameter(cyt, initial=0)
reaction_gaba_flux = rxd.Rate(gaba, gaba_flux_rate)
#Set the reactions between the species
reaction000 = rxd.Reaction(gaba, gabaOut, ks[0])
reaction001 = rxd.Reaction(gaba + GABABR, gabaGABABR, ks[1], ks[2])
reaction002 = rxd.Reaction(gabaGABABR + Gi, gabaGABABRGi, ks[3], ks[4])
reaction003 = rxd.Reaction(gabaGABABRGi, gabaGABABRGibg + GiaGTP, ks[5])
reaction004 = rxd.Reaction(gabaGABABRGibg, gabaGABABR + Gibg, ks[6])
reaction005 = rxd.Reaction(GABABR + Gi, GABABRGi, ks[7], ks[8])
reaction006 = rxd.Reaction(gaba + GABABRGi, gabaGABABRGi, ks[9], ks[10])
reaction007 = rxd.Reaction(GiaGTP + RGS, GiaGTPRGS, ks[11], ks[12])
reaction008 = rxd.Reaction(GiaGTPRGS, GiaGDP + RGS, ks[13])
reaction009 = rxd.Reaction(GiaGDP + Gibg, Gi, ks[14])
reaction010 = rxd.Reaction(GIRK + Gibg, GIRKGibg, ks[15], ks[16])
reaction011 = rxd.Reaction(GIRKGibg + Gibg, GIRKGibg2, ks[17], ks[18])
reaction012 = rxd.Reaction(GIRKGibg2 + Gibg, GIRKGibg3, ks[19], ks[20])
reaction013 = rxd.Reaction(GIRKGibg3 + Gibg, GIRKGibg4, ks[21], ks[22])
reaction014 = rxd.Reaction(VGCC + Gibg, VGCCGibg, ks[23], ks[24])
reaction015 = rxd.Reaction(4*ca + Ffactor, Ffactorca4, ks[25], ks[26]) #In Destexhe 1994 this reaction doesn't affect Ca concentration, here it does
reaction016 = rxd.Reaction(Ffactorca4 + Ves, Ffactorca4Ves, ks[27], ks[28])
reaction017 = rxd.Reaction(Ffactorca4Ves, 10000*NT + Ffactor + Ves + 4*ca, 10000*ks[29]*Ffactorca4Ves, 0, custom_dynamics=True) #Note: this is a bit different than in Destexhe 1994. Here the amount of Vesicles is dynamical, in Destexhe Ves was constant. Here the used vesicle and ca returns to baseline
reaction018 = rxd.Reaction(NT, NT*0, ks[30])
reaction019 = rxd.Rate(ca, -(ca - minCa)/decay)
#Recordings:
vec_t = h.Vector()
vec_t.record(h._ref_t)
vecs = []
for ispec in range(0,len(species)):
vecs.append(h.Vector())
print('First recording...')
timenow = time.time()
vecs[0].record(gaba.nodes(h.spine)(0.5)[0]._ref_concentration)
print('First recording set in '+str(time.time()-timenow)+' seconds')
vecs[1].record(gabaOut.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[2].record(GABABR.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[3].record(gabaGABABR.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[4].record(Gi.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[5].record(GABABRGi.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[6].record(gabaGABABRGi.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[7].record(gabaGABABRGibg.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[8].record(GiaGTP.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[9].record(GiaGDP.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[10].record(RGS.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[11].record(GiaGTPRGS.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[12].record(Gibg.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[13].record(GIRK.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[14].record(VGCC.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[15].record(GIRKGibg.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[16].record(GIRKGibg2.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[17].record(GIRKGibg3.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[18].record(GIRKGibg4.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[19].record(VGCCGibg.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[20].record(ca.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[21].record(Ffactor.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[22].record(Ffactorca4.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[23].record(Ves.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[24].record(Ffactorca4Ves.nodes(h.spine)(0.5)[0]._ref_concentration)
vecs[25].record(NT.nodes(h.spine)(0.5)[0]._ref_concentration)
print('Biochem recordings OK')
axon = h.axon
spine = h.spine
cyt = rxd.Region([spine], name='cyt', nrn_region='i')
#Set the stimulations and the recordings
h("""
objref st1s
st1s = new List()
for (i=0;i<"""+str(Nstims-1)+""";i+=1) {
axon[0] st1s.append(new IClamp(0.5))
st1s.o[i].del = 1000 + """+str(stim_T)+"""*(i+1)
st1s.o[i].dur = 0.5
st1s.o[i].amp = """+str(stimAmp)+"""
}
objref vrec, tvec, carec, icarec, Trec, VArec, FArec
vrec = new Vector()
carec = new Vector()
icarec = new Vector()
Trec = new Vector()
VArec = new Vector()
FArec = new Vector()
tvec = new Vector()
spine cvode.record(&v(0.5),vrec,tvec)
spine cvode.record(&cai(0.5),carec,tvec)
spine cvode.record(&ica_CaNtypeWillingReluctant(0.5),icarec,tvec)
//spine cvode.record(&T_rel(0.5),Trec,tvec)
//spine cvode.record(&VA_rel(0.5),VArec,tvec)
//spine cvode.record(&FA_rel(0.5),FArec,tvec)
//spine icarec.record(&ica_CaNtypeWillingReluctant(0.5),0.5)
cai0_ca_ion = """+str(minCa)+"""
cao0_ca_ion = 2
spine area_spine = area(0.5)
celsius = 33
""")
print('Electrophys recordings OK')
cols = mytools.colorsredtolila(7)
mycol = cols[0]
A_spine = h.area_spine*1e-8 #multiplied by 1e-8 for conversion um2 -> cm2
iVGCC = [i for i in range(0,len(species)) if species[i] == 'VGCC'][0] #14
N_VGCC = initvalues[iVGCC]*1e-3*my_volume*6.022e23 #mM -> M -> mol -> number of molecules
gCaNtypebar_tot = N_VGCC*gCaNtypebar_single #uS
gCaNtypebar_per_area = gCaNtypebar_tot/A_spine*1e-6 #uS -> uS/cm2 -> S/cm2
h('spine gcalbar_CaNtypeWillingReluctant = '+str(gCaNtypebar_per_area))
timenow = time.time()
h.init()
#Set the events (GABA inputs)
def set_param(param, val):
param.nodes.value = val
h.cvode.re_init()
def add_param(param, val):
print("Old param.nodes.value = "+str(param.nodes.value)+", t="+str(h.t))
param.nodes.value = [q + val for q in param.nodes.value]
print("New param.nodes.value = "+str(param.nodes.value)+", t="+str(h.t))
h.cvode.re_init()
h.cvode.event(1000, lambda: add_param(gaba_flux_rate, gaba_input_flux/6.022e23/my_volume*1e3))
h.cvode.event(1000+gaba_input_dur, lambda: add_param(gaba_flux_rate, -gaba_input_flux/6.022e23/my_volume*1e3))
h.continuerun(tstop)
print("Simulation gCaNtypebar_per_area="+str(gCaNtypebar_per_area)+" done in "+str(time.time()-timenow)+" seconds")
#Plot and save the results
times = array(h.tvec)
Vrec= array(h.vrec)
Carec= array(h.carec)
iCarec= array(h.icarec)
VGCCrec = array(vecs[14])
VGCCGibgrec = array(vecs[19])
FArec= array(vecs[22])
VArec= array(vecs[24])
Trec= array(vecs[25])
f,axarr = subplots(3,3)
for iax in range(0,3):
for iay in range(0,3):
axarr[iay,iax].tick_params(axis='both', which='major', labelsize=4)
axarr[0,0].set_title('$V_m$, gnabar='+str(gnabar),fontsize=5)
axarr[1,0].set_title('Ca2+ current (mA/cm2), gkbar='+str(gkbar),fontsize=5)
axarr[2,0].set_title('[Ca2+], gl='+str(gl),fontsize=5)
axarr[0,1].set_title('[FA] (mM)',fontsize=5)
axarr[1,1].set_title('[VA] (mM)',fontsize=5)
axarr[2,1].set_title('[neurotrans.] (mM)',fontsize=5)
axarr[0,2].set_title('[gaba] (mM)',fontsize=5)
axarr[1,2].set_title('[Gibg] (mM)',fontsize=5)
axarr[2,2].set_title('[VGCCGibg] (mM)',fontsize=5)
axarr[0,0].plot(times,Vrec,color=mycol,lw=0.5)
axarr[1,0].plot(times,iCarec,color=mycol,lw=0.5)
axarr[2,0].plot(times,Carec,color=mycol,lw=0.5)
it = argmax(Carec)
it_minus_10ms = argmin(abs(times-times[it]+10))
it_plus_50ms = argmin(abs(times-times[it]-50))
it_max50perc_before = argmin([abs(Carec[i]-0.5*max(Carec))+1e5*(i<it_minus_10ms or i>=it) for i in range(0,len(times))])
it_max50perc_after = argmin([abs(Carec[i]-0.5*max(Carec))+1e5*(i>=it_plus_50ms or i<it) for i in range(0,len(times))])
print("Max Ca at t="+str(times[argmax(Carec)])+", ca transient width = "+str(times[it_max50perc_after]-times[it_max50perc_before])) #t_max50perc_before="+str(times[it_max50perc_before])+", t_max50perc_after="+str(times[it_max50perc_after])+", max = "+str(max(Carec))+", after="+str(Carec[it_max50perc_after]))
axarr[0,1].plot(times,FArec,color=mycol,lw=0.5)
axarr[1,1].plot(times,VArec,color=mycol,lw=0.5)
axarr[2,1].plot(times,Trec,color=mycol,lw=0.5)
iTpeak = argmax(Trec)
iThalfpeak = argmin(abs(Trec-0.5*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iThalfperpeak = argmin(abs(Trec-0.5/exp(1)*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT2 = argmin(abs(Trec-0.75*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT3 = argmin(abs(Trec-0.75/exp(1)*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT4 = argmin(abs(Trec-0.25*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT5 = argmin(abs(Trec-0.25/exp(1)*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT6 = argmin(abs(Trec-0.2*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT7 = argmin(abs(Trec-0.2/exp(1)*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT8 = argmin(abs(Trec-0.15*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT9 = argmin(abs(Trec-0.15/exp(1)*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT10 = argmin(abs(Trec-0.1*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
iT11 = argmin(abs(Trec-0.1/exp(1)*max(Trec)) - 1000*(times > times[iTpeak])*(times < times[iTpeak]+500)) #allow time points from the peak until 500 ms after the peak, other time points are given 1000 extra error
print("NT max "+str(max(Trec))+", decay constant = "+str(times[iThalfperpeak] - times[iThalfpeak])+" or "+str(times[iT3]-times[iT2])+" or "+str(times[iT5]-times[iT4]))
axarr[0,2].plot(times,vecs[0],color=mycol,lw=0.5)
axarr[1,2].plot(times,vecs[12],color=mycol,lw=0.5)
axarr[2,2].plot(times,vecs[19],color=mycol,lw=0.5)
axarr[2,2].plot(times,vecs[14])
for i in range(0,3):
axarr[i,0].set_xlim([900,1000+stim_T+1200])
axarr[i,1].set_xlim([900,1000+stim_T+1200])
axarr[i,2].set_xlim([900,1000+stim_T+1200])
addition = '_EE' if isyntype == 1 else ('_IE' if isyntype == 2 else '_wrongRGS')
if len(blockedStr) > 0 and blockedStr != 'None':
addition = addition+'_'+blockedStr
axarr[2,1].set_ylim([0,5.5])
if doDraw:
f.savefig("presynGABA_rxdrel6_no1stglu_recVGCC"+addition+"_gnabar"+str(gnabar)+"_gkbar"+str(gkbar)+"_gl"+str(gl)+"_el"+str(el)+"_decay"+str(decay)+'_minCa'+str(minCa)+"_"+params_str+"_Ves"+str(Vesmax)+"_Fmax"+str(Fmax)+"_"+str(stim_T)+"ms_N"+str(Nstims)+"_len"+str(axonPasLen)+","+str(axonActLen)+'_amp'+str(stimAmp)+'_kcoeff'+str(kcoeff_f)+"_"+str(kcoeff_b)+"_"+str(kcoeff_2f)+"_"+str(kcoeff_2b)+"_"+str(kcoeff_3)+"_"+str(kcoeff_4)+".pdf")
filename = "presynGABA_rxdrel6_no1stglu_recVGCC"+addition+"_gnabar"+str(gnabar)+"_gkbar"+str(gkbar)+"_gl"+str(gl)+"_el"+str(el)+"_decay"+str(decay)+'_minCa'+str(minCa)+"_"+params_str+"_Ves"+str(Vesmax)+"_Fmax"+str(Fmax)+"_"+str(stim_T)+"ms_N"+str(Nstims)+"_len"+str(axonPasLen)+","+str(axonActLen)+'_amp'+str(stimAmp)+'_kcoeff'+str(kcoeff_f)+"_"+str(kcoeff_b)+"_"+str(kcoeff_2f)+"_"+str(kcoeff_2b)+"_"+str(kcoeff_3)+"_"+str(kcoeff_4)+'.mat'
if fileOut != 'None':
filename = fileOut+'.mat'
if not doInterpolate:
scipy.io.savemat(filename,{'times': times, 'NT': Trec,'Vrec': Vrec, 'Carec': Carec, 'VGCCGibgrec': VGCCGibgrec, 'VGCCrec': VGCCrec})
else:
ts = [1.0*i for i in range(0,int(tstop))]
scipy.io.savemat(filename,{'times': ts, 'NT': mytools.interpolate_extrapolate_constant(times,Trec,ts),'Vrec': mytools.interpolate_extrapolate_constant(times,Vrec,ts), 'Carec': mytools.interpolate_extrapolate_constant(times,Carec,ts), 'VGCCrec': mytools.interpolate_extrapolate_constant(times,VGCCrec,ts), 'VGCCGibgrec': mytools.interpolate_extrapolate_constant(times,VGCCGibgrec,ts), 'FArec': mytools.interpolate_extrapolate_constant(times,FArec,ts), 'VArec': mytools.interpolate_extrapolate_constant(times,VArec,ts)})