// Octopus cell, stick representation
// From Spencer et al., Front. Comput. Neurosci., 22 October 2012 | https://doi.org/10.3389/fncom.2012.00083
// Figure 1.
load_file("nrngui.hoc")
load_proc("nrnmainmenu")

objref soma
soma = new SectionList()
objref syn
objref ncl, weight
objref apc
strdef name
objref primarydendrite
primarydendrite = new SectionList()
objref hillock
hillock = new SectionList()
objref unmyelinatedaxon
unmyelinatedaxon = new SectionList()

tstop = 1000

RM= 5000   // ohm cm2
ra=150  // ohm cm
c_m=0.9  //muF cm2
resting_Vm=-65
totcap =12 // total cap in pF for cell 
somaarea = totcap *1E-6 / 1 // pf -> uF,assumes 1 uF/cm2; result is in cm2 
lstd = 1E4*sqrt(somaarea/PI) // convert from cm to um 

/*
gna=0.083333333                  //gna=0.07958
gh=0.005                    //gh=0.0005
gka=0.00477                  //gka= 0.00477
gkht=0.1875                //gkht=0.01592
gklt=0.75                //gklt=0.01592
g_leak=0.0016666667
*/

create sections[7]

access sections[0]
soma.append()
sections[0] {
    pt3dadd(0., 0., 0., 25.)
    pt3dadd(0., 25., 0., 25.)
}

// axon hillock, 30 microns long, untapered
access sections[1]
hillock.append()
connect sections[1](0), sections[0](0)
sections[1] {
    pt3dadd(0., 0., 0., 3)
    pt3dadd(0., -30., 0., 3)
}

// initial segment, unmyelinated, 10 microns long, not tapering
access sections[2]
unmyelinatedaxon.append()
connect sections[2](0), sections[1](1)
sections[2] {
    pt3dadd(0., -30., 0., 3.0)
    pt3dadd(0., -32., 0., 3.0)
}

// angle dendrites out a bit from original model for visibility, so 
// length for 20 u end deviation is 279.285 um, end pos 304.285

access sections[3]
primarydendrite.append()
connect sections[3](0), sections[0](1)
sections[3] {
    pt3dadd(0., 25., 0., 3)
    pt3dadd(20., 304.285, 0., 3)
}

access sections[4]
primarydendrite.append()
connect sections[4](0), sections[0](1)
sections [4] {
    pt3dadd(0., 25., 0., 3)
    pt3dadd(0., 304.285, 20., 3)
}

access sections[5]
primarydendrite.append()
connect sections[5](0), sections[0](1)
sections [5] {
    pt3dadd(0., 25., 0., 3)
    pt3dadd(-20., 304.285, 0., 3)
}

access sections[6]
primarydendrite.append()
connect sections[6](0), sections[0](1)
sections [6] {
    pt3dadd(0., 25., 0., 3)
    pt3dadd(0., 304.285, -20., 3)
}

for i=0,6 {
sections[i] {
    insert pas
    Ra=ra
    cm=c_m
    e_pas=resting_Vm
    g_pas=1/RM
}
}

// convert from nanosiemens to mho/cm2.
func nstomho() {
    return (1E-9*$1/somaarea)
}


sections[1] {
    insert jsrnaf ena_jsrnaf = 50
    insert klt ek_lt = -70
    insert kht ek_kht = -70
    insert na ena_na = 50
    insert ka ek_ka=-70
    insert ih eh_ih=-43
    insert hcno eh_hcno = -43 
    insert leak g_leak=1/10000 erev_leak = -65

    gnabar_jsrnaf = 10*nstomho(1000)
    gkhtbar_kht = 2*nstomho(150)
    gkltbar_klt = 2*nstomho(600)
    gkabar_ka = 150*nstomho(200)
    ghbar_ih = nstomho(0)
    gbar_hcno = 5*nstomho(40)
    g_leak = 10*nstomho(2)

    print " " , gnabar_jsrnaf
    print " " , gkltbar_klt
    print " " , gkhtbar_kht
    print " " , g_leak
    print " " , gbar_hcno
}

sections[2] {
    insert jsrnaf ena_jsrnaf = 50
    insert klt ek_lt = -70
    insert kht ek_kht = -70
    insert na ena_na = 50
    insert ka ek_ka=-70
    insert ih eh_ih=-43
    insert hcno eh_hcno = -43 gbar_hcno = 0
    insert leak g_leak=1/10000 erev_leak = -65


    gnabar_jsrnaf = 10*nstomho(1000)
    gkhtbar_kht = 2*nstomho(150)
    gkltbar_klt = 2*nstomho(600)
    gkabar_ka = 150*nstomho(200)
    ghbar_ih = nstomho(0)
    gbar_hcno = 2.5*nstomho(40)
    g_leak = 10*nstomho(2)
}

sections[0] {
    insert jsrnaf ena_jsrnaf = 50
    insert klt ek_lt = -70
    insert kht ek_kht = -70
    insert na ena_na = 50
    insert ka ek_ka=-70
    insert ih eh_ih=-43
    insert hcno eh_hcno = -43 
    insert leak g_leak=1/10000 erev_leak = -65

    gnabar_jsrnaf = 10*nstomho(1000)
    gkhtbar_kht = 2*nstomho(150)
    gkltbar_klt = 2*nstomho(600)
    gkabar_ka = 150*nstomho(200)
    ghbar_ih = nstomho(0)
    gbar_hcno = 2.5*nstomho(40)
    g_leak = 10*nstomho(2)
}


for i=3,6 {
sections[i] {
    insert klt ek_lt = -70
    insert kht ek_kht = -70
    insert ka ek_ka=-70
    insert ih eh_ih=-43
    insert hcno eh_hcno = -43 
    insert leak g_leak=1/10000 erev_leak = -65

    gkhtbar_kht = 2*nstomho(150)
    gkltbar_klt = 2*nstomho(600)
    gkabar_ka = 150*nstomho(200)
    ghbar_ih = nstomho(0)
    gbar_hcno = 2.5*nstomho(40)
    g_leak = 10*nstomho(2)
}
}



objref syn[500], psyn[500], stim, stim1
nsyn=1
gmax=0 //uS (Minimum:0.5nS and Maximum 35nS)
/*
sections[6] { 
    syn[0]= new syn_g_Octopus(1)
    syn[0].onset=100
    syn[0].gmax=gmax
}
*/

amp=10
dur=150

amp1=4
dur1=20

sections [0] stim = new IClamp(0.5)

stim.del = 400
stim.dur = dur
stim.amp = amp

/*
sections [0] stim1 = new IClamp(0.5)

stim1.del = 420
stim1.dur = dur1
stim1.amp = amp1
*/


proc update_CM() {

    for j=0, 6 {
        sections[i] {
            cm=c_m
        }
    }    

}
proc update_ra() {

     for j=0, 6 {
        sections[i] {
             Ra=ra
        }
     }   
}
proc update_RM() {

      for j=0, 6 {
        sections[i] {
        g_pas=1/RM
        }
    }
}

proc update_RestingVm() {

    for j=0, 6 {
        sections[i] {
        e_pas=resting_Vm
        }
    }    
    
}

xpanel("Parameters")

xvalue("Rm","RM",1,"update_RM()",1,1)
xvalue("Ra","ra",1,"update_ra()",1,1)
xvalue("Cm","c_m",1,"update_CM()",1,1)
xvalue("Resting_Vm","resting_Vm",1,"update_RestingVm()",1,1)

xpanel()

objref g[20]
ngraph=0

proc addgraph() { local ii  // define subroutine to add a new graph
                // addgraph("variable", minvalue, maxvalue)
    ngraph = ngraph+1
    ii = ngraph-1
    g[ii] = new Graph()
    g[ii].size(0,tstop,$2,$3)
    g[ii].xaxis()
    g[ii].yaxis()
    g[ii].addvar($s1,1,0)
    g[ii].save_name("graphList[0].")
    graphList[0].append(g[ii])
}


objref stfunc,shape
shape=new Shape(0)
shape.view(-792.112, -154.611, 1363.32, 1185.94, 5, 340, 300, 260)
shape.show(0)
proc make_shape_plot(){//DRAWS THE POINTS ON THE CELL
    shape.point_mark_remove()

    for i=0,nsyn{
        shape.point_mark(syn[i], 2, 4, 4)   
    }

}///END SHAPE

addgraph("sections[0].v(0.5)",-100,100)
addgraph("sections[6].v(0.5)",-100,100)
//addgraph("syn[0].inmda",-1,1)
//addgraph("dend.v(0.5)",-100,100)
//addgraph("soma.i_pas(0.5)",-100,100)
//g[0].addexpr("soma.v(0)",2,1)
//g[1].addexpr("syn[1].gampa",-1,1)
//g[2].addexpr("syn[1].iampa",-1,1)
make_shape_plot()

nrncontrolmenu()