spikegen_GPeA.hoc

{load_file("nrngui.hoc")}
{load_file("stdlib.hoc")}

//DEBUG_MSG = 1

{load_file("ParaCellSwc.hoc")}
{load_file("areaList.hoc")}
{load_file("exSectionList.hoc")}
/*
USEAGE
./x86_64/special -c "strdef FILENAME" -c "FILENAME=voltage00.txt" -c "SPIKEGEN_POS=70" -c "SPIKEGEN_SIZE=10" -c "SPIKEGEN_GK=0.36" -c "SPIKEGEN_GNA=1.2" -c "SPIKEGEN_GCA=0.02" -c "DEBUG_MSG=1" spikegen_model.hoc

*/

if(name_declared("MODEL")!=5 ){ execute("MODEL=6") }
if(name_declared("FILE_SUFFIX")!=5){ execute("FILE_SUFFIX=9999")}
if(name_declared("SPIKEGEN_POS")!=5){ execute("SPIKEGEN_POS=10")}
if(name_declared("SPIKEGEN_SIZE")!=5){ execute("SPIKEGEN_SIZE=10") }
if(name_declared("GKDRBAR")!=5){ execute("GKDRBAR=0.7") }
if(name_declared("GNABAR")!=5){ execute("GNABAR=0.5") }
if(name_declared("GCATBAR")!=5){ execute("GCATBAR=0.0015") }
if(name_declared("GKCABAR")!=5){ execute("GKCABAR=0.15") }
if(name_declared("KCA")!=5){ execute("KCA=2.4") }
if(name_declared("KNA")!=5){ execute("KNA=5.0") }
if(name_declared("KK")!=5){ execute("KK=KNA") }
if(name_declared("CM")!=5){ execute("CM=1") } //
if(name_declared("GL")!=5){ execute("GL=0.00035") } //
if(name_declared("EL")!=5){ execute("EL=-60") } //
if(name_declared("CAGAIN")!=5){ execute("CAGAIN=0.16") }
if(name_declared("DEBUG_MSG")!=5){ execute("DEBUG_MSG=1") }
if(!name_declared("FILENAME")){ execute("strdef FILENAME FILENAME=\"targetdata11.txt\"") }

/*
if(SET_VALUE==1){
    strdef FILENAME
    FILENAME="output/voltage.txt"

    SPIKEGEN_POS=70
    SPIKEGEN_SIZE=10
    SPIKEGEN_GK=0.036*10
    SPIKEGEN_GNA=0.12*10
    SPIKEGEN_GCA=0.002*10
    DEBUG_MSG=1
}
*/
//if(name_declared("FILENAME")!=4 || strcmp(FILENAME,"")){ FILENAME="output/voltage.txt" printf("test") }
//printf("filename=%s\n", FILENAME)
//printf("dec=%d\n", name_declared("FILENAME"))
//printf("dec=%d\n", name_declared("FILENAME"))

NCELL = 1
STOPTIME = 200
//STIM_POINT = 1106
//MEASURE_POINT1 = 1100
//MEASURE_POINT2 = 1099
STIM_POINT = 9
MEASURE_POINT1 = 1208
MEASURE_POINT2 = 948


/*
objref pc
pc = new ParallelContext()
*/

/*
SPIKEGEN_POS = 70
SPIKEGEN_SIZE = 10
SPIKEGEN_GK = 0.036*10
SPIKEGEN_GNA = 0.12*10
SPIKEGEN_GCA = 0.002*10
*/

//{printf("pos=%d, size=%d, G_K=%f, G_Na=%f, G_Ca=%f\n", SPIKEGEN_POS, SPIKEGEN_SIZE, SPIKEGEN_GK, SPIKEGEN_GNA, SPIKEGEN_GCA)}

//objref cells, nclist
objref nclist
objref g_ncstim_list
objref stim, ncstim, stim_vec, vec_t
objref g_spike_time_vec, g_spike_id_vec
objref g_volt_time_vec, g_volt_stim_vec, g_volt_obs_vec1, g_volt_obs_vec2
objref g_cur_stim_vec
objref spikegenlist, noaxonlist

proc makeCells(){local i, axonnum localobj cell, nc, nil, tmplist
    /*
    cells = new List()
    for(i=pc.id; i<$1; i +=pc.nhost){
	cell = new CellSwc("./swc/PN/toroid/050622_4_sn.swc")

	cells.append(cell)
	pc.set_gid2node(i, pc.id)
	nc = cell.makeNetCon(0, nil, 0, 0, 1.0)
	pc.cell(i, nc)
    }
    */
    p_mkcells("neuronlist.txt",0)
    cellbodylist = areaList(39,38,0) // gyaku?
    axonlist = areaList(99,100,0)
    dendlist = areaList(97,98,0)
    noaxonlist = new SectionList()

    cells.o(0).Dend[0] noaxonlist.wholetree()
    noaxonlist.remove(axonlist)

    spikegenlist = new SectionList()
    axonnum = countSL(axonlist)
    tmplist = sl2list(axonlist)
    i = SPIKEGEN_POS
    while(i < axonnum && i < SPIKEGEN_POS+SPIKEGEN_SIZE) {
	tmplist.o(i).sec spikegenlist.append()
	i+=1
    }
    /*
    cells.o(0).Dend[SPIKEGEN_POS] stmp1 = secname()
    stmp2 = ""
    tmplist = new List()
    forsec axonlist tmplist.append(new SectionRef())
    i = 0
    while( strcmp(stmp1,stmp2) != 0 ) {
	tmplist.o(i).sec stmp2 = secname()
	i+=1
    }
    for(j=i-1; j<i+SPIKEGEN_SIZE-1; j+=1) {
	if( j < tmplist.count() ) {
	    tmplist.o(j).sec spikegenlist.append()
	}
    }
    */
    forsec noaxonlist {
	insert pas
	e_pas = EL
	g_pas = GL
	cm = CM
    }
    /*
    forall {
	insert pas
	e_pas = EL
	g_pas = GL
    }
    */
    forsec axonlist {
	insert GPeA_k
	gnabar_GPeA_k = GNABAR*0.1	 
	gkdrbar_GPeA_k = GKDRBAR*0.1
	gl_GPeA_k = GL
	el_GPeA_k = EL
	gcatbar_GPeA_k = GCATBAR*0.1
	gkcabar_GPeA_k = GKCABAR*0.1
	kca_GPeA_k = KCA
	kna_GPeA_k = KNA
	kk_GPeA_k = KNA
	caGain_GPeA_k = CAGAIN
	vol_GPeA_k = PI*diam*diam*L/4 * 1e-15
	cm = CM
    }
    
    forsec spikegenlist {
	gnabar_GPeA_k = GNABAR
	gkdrbar_GPeA_k = GKDRBAR
	gcatbar_GPeA_k = GCATBAR
	gkcabar_GPeA_k = GKCABAR
    }
	

    /*
    for(i=0; i<102; i=i+1){
	CellSwc[0].Dend[i]{
            insert hhext2
	    insert CaT
	    gkbar_hhext2 = 0.036
	    gnabar_hhext2 = 0.12
	    gmax_CaT = 0.002
	}
    }
    for(i=102; i<1112; i=i+1){
	CellSwc[0].Dend[i]{
	    insert pas
	}
    }
    for(i=SPIKEGEN_POS; (i<SPIKEGEN_POS+SPIKEGEN_SIZE && i<102); i=i+1){
	CellSwc[0].Dend[i]{
 	    gkbar_hhext2 = SPIKEGEN_GK
	    gnabar_hhext2 = SPIKEGEN_GNA
	    gmax_CaT = SPIKEGEN_GCA
	}
    }
    */
}


proc setStim() { localobj cell, synStim
    if(!pc.gid_exists(0)) { return }
    stim = new NetStim()
    stim.number = 1
    stim.start = 0.5
    cell = pc.gid2cell(0)
    synStim = cell.makeSyn(STIM_POINT, 2, 5, 0)

    ncstim = new NetCon(stim, synStim)
    ncstim.delay = 0.1
    ncstim.weight = 0.01
}

proc setIClamp() { local span, amp
    if(!pc.gid_exists(0)) { return }
    
    amp = 5
    span = STOPTIME/dt
    cells.o(0).Dend[STIM_POINT] stim = new IClamp(0.5)
    vec_t = new Vector(span)
    stim_vec = new Vector(span)
    for(i=0; i<span; i=i+1){ vec_t.x(i) = dt*i }
    //for(i=0; i<200/dt; i=i+1){ stim_vec.x(i) = 10*i/(200/dt) }
    stim.del = 0
    stim.dur = 1e9
    if(MODEL==4 || MODEL==7){
	for(i=0; i<200/dt; i=i+1){ stim_vec.x(i) = amp*i/(200/dt) }
    }else if(MODEL==5 || MODEL==8){
	for(i=0; i<200/dt; i=i+1){ stim_vec.x(i) = amp*0.1*sin(5*PI*i/(100/dt)) + amp/2*i/(200/dt) }
	//for(i=0; i<200/dt; i=i+1){ stim_vec.x(i) = 1.0*sin(5*PI*i/(100/dt)) + 5*i/(200/dt) }
    }else if(MODEL==6 || MODEL==9){
	for(i=0; i<200/dt; i=i+1){ stim_vec.x(i) = amp*0.1*sin(10*PI*i/(100/dt)) + amp*0.5*i/(200/dt) }
	//for(i=0; i<200/dt; i=i+1){ stim_vec.x(i) = 1.0*sin(10*PI*i/(100/dt)) + 5*i/(200/dt) }
    }
    
    if(MODEL==10){
	stim.amp = amp
    }else{
	stim_vec.play(&stim.amp, vec_t, 1)
    }

    //stim.dur = 200
    //stim.amp = 10
}

proc setIClamp2() { local span localobj f
    if(!pc.gid_exists(0)) { return }
    strdef tmp
    
    dt = 0.05
    steps_per_ms = 20
    //span = STOPTIME/dt
    cells.o(0).Dend[STIM_POINT] stim = new IClamp(0.5)

    f = new File()
    { f.ropen(FILENAME) }
    for (i=0; i<2; i=i+1) { f.gets(tmp) }
    span = f.scanvar()
    { f.scanvar() }
    
    vec_t = new Vector(span)
    stim_vec = new Vector(span)
    for (i=0; i<span; i=i+1) {
	vec_t.x(i) = f.scanvar()
	stim_vec.x(i) = f.scanvar()
	{f.scanvar()}
    }
    f.close()
    
    stim.del = 0
    stim.dur = 1e9
    {stim_vec.play(&stim.amp, vec_t, 1)}
    
    STOPTIME = int(span * vec_t.x(1))
}

proc setSpikeRecord() { local i localobj nc0, nc1, nc2, nil
    g_spike_time_vec = new Vector()
    g_spike_id_vec = new Vector()
    g_ncstim_list = new List()
    
    nc0 = cells.object(i).makeNetCon(STIM_POINT, nil, 0, 0, 1.0)
    pc.set_gid2node(100+1, pc.id)
    pc.cell(100+1, nc0)
    nc0.record(g_spike_time_vec, g_spike_id_vec, nc0.srcgid)
    g_ncstim_list.append(nc0)
    
    nc1 = cells.object(i).makeNetCon(MEASURE_POINT1, nil, 0, 0, 1.0)
    pc.set_gid2node(100+2, pc.id)
    pc.cell(100+2, nc1)
    nc1.record(g_spike_time_vec, g_spike_id_vec, nc1.srcgid)
    g_ncstim_list.append(nc1)

    nc2 = cells.object(i).makeNetCon(MEASURE_POINT2, nil, 0, 0, 1.0)
    pc.set_gid2node(100+3, pc.id)
    pc.cell(100+3, nc2)
    nc2.record(g_spike_time_vec, g_spike_id_vec, nc2.srcgid)
    g_ncstim_list.append(nc2)
}

proc setVoltageRecord() {
    g_volt_time_vec = new Vector()
    g_volt_stim_vec = new Vector()
    g_volt_obs_vec1 = new Vector()
    g_volt_obs_vec2 = new Vector()
    
    g_volt_time_vec.record(&t)
    //g_volt_stim_vec.record(&v(0.5))
    g_volt_stim_vec.record( &cells.o(0).Dend[STIM_POINT].v(0.5) )
    g_volt_obs_vec1.record( &cells.o(0).Dend[MEASURE_POINT1].v(0.5) )
    g_volt_obs_vec2.record( &cells.o(0).Dend[MEASURE_POINT2].v(0.5) )
}

proc setCurrentRecord() {
    g_cur_stim_vec = new Vector()
    g_cur_stim_vec.record(&stim.i)
}

proc printVoltage(){ localobj fobj, mt
    pc.barrier()
    
    mt = new Matrix()
    printf("#  t (msec)\tIClamp\tVobs1 (%d)\n", STIM_POINT)
    mt.resize(g_volt_time_vec.size(), 3)
    mt.setcol(0, g_volt_time_vec)
    mt.setcol(1, g_cur_stim_vec)
    if(MODEL==4 || MODEL==5 || MODEL==6){
	mt.setcol(2, g_volt_obs_vec1)
    }else if(MODEL==7 || MODEL==8 || MODEL==9 || MODEL==10){
	mt.setcol(2, g_volt_obs_vec2)
    }
    mt.printf("%0.5f ")
}

proc writeVoltage() { localobj fobj, mt
    pc.barrier()
    strdef filename
    
    fobj = new File()
    mt = new Matrix()	
    //filename = new String()
    sprint(filename, "output/voltage%010d.txt", FILE_SUFFIX)
    fobj.wopen(filename)
    fobj.printf("  t (msec)\tVobs1 (%d)\n", STIM_POINT, MEASURE_POINT1, MEASURE_POINT2)
    //fobj.printf("  t (msec)\tVstim (%d)\tVobs1 (%d)\tVobs2 (%d)\n", STIM_POINT, MEASURE_POINT1, MEASURE_POINT2)
    mt.resize(g_volt_time_vec.size(), 3)
    mt.setcol(0, g_volt_time_vec)
    mt.setcol(1, g_cur_stim_vec)
    //mt.setcol(2, g_volt_stim_vec)
    if(MODEL==4 || MODEL==5 || MODEL==6){
	mt.setcol(2, g_volt_obs_vec1)
    }else if(MODEL==7 || MODEL==8 || MODEL==9 || MODEL==10){
	mt.setcol(2, g_volt_obs_vec2)
    }
    //mt.setcol(2, g_volt_obs_vec1)
    //mt.setcol(4, g_volt_obs_vec2)
    mt.fprint(0, fobj, "  %0.5f")
    fobj.close()
}

proc writeSpike(){ localobj fobj, mt
    pc.barrier()
    strdef filename
    
    fobj = new File()
    mt = new Matrix()	
    //	filename = new String()
    sprint(filename, "output/output%010d.txt", pc.id)
    fobj.wopen(filename)
    fobj.printf("time [ms]\t cell_id\n")
    mt.resize(g_spike_time_vec.size(), 2)
    mt.setcol(0, g_spike_time_vec)
    mt.setcol(1, g_spike_id_vec)
    mt.fprint(0, fobj, "\t%g\t")
    fobj.close()
}

proc printSpike() { local i, rank
    pc.barrier()
    if(pc.id==0) printf("\ntime\t cell\n")
    for rank=0, pc.nhost-1 {
	if(rank==pc.id){
	    for i=0, g_spike_time_vec.size-1 {
		printf("%g\t %d\n", g_spike_time_vec.x[i], g_spike_id_vec.x[i])
	    }
	}
	pc.barrier()
    }
}


proc main() {local i, start_tiem,modelfin_time, stop_time, start_waittime
    start_time = pc.time

    makeCells(1)
    //connectcells()

    //setStim()
    //setIClamp2()
    setIClamp()
    setSpikeRecord()
    setVoltageRecord()
    setCurrentRecord()
    tstop = STOPTIME
    pc.set_maxstep(10)
    stdinit()
    
    pc.barrier()
    modelfin_time = pc.time
    start_waittime = pc.wait_time
    
    //forall{ psection() }
    //topology()
    
    pc.psolve(tstop)
    pc.barrier()
    stop_time = pc.time
    if(pc.id ==0 && DEBUG_MSG==1){
	printf("[Process=%d, Thread=%d]\n", pc.nhost, pc.nthread)
        printf(" * NCELL=%d, NSTEP=%d (STOPTIME=%f / dt=%f)\n", NCELL, STOPTIME/dt, STOPTIME, dt)
	printf(" * step : %f sec, wait : %f sec, send : %f sec\n", pc.step_time(), pc.wait_time()-start_waittime, pc.send_time)
	printf(" * modeling time : %f sec\n", modelfin_time - start_time)
        printf(" * core time : %f sec\n", stop_time-modelfin_time)
	
    }
        
    //printSpike()
    //writeRecord()
    //writeVoltage()
    printVoltage()
    {pc.runworker()}
    {pc.done()}

}

main()
{quit()}