load_file("nrngui.hoc")

/* simulation of neuronal model with varying complexity (set nr 'branch_levels').
A constant apical trunk diameter is used and the rall diameter rule is applied

license: 
author: Luuk van der Velden, University of Amsterdam
paper: Altered dendritic complexity affects firing properties of cortical layer 2/3 pyramidal neurons in mice lacking the 5-HT3A receptor
journal of neurophysiology, 2012

no warranties :)

run as: 'nrngui filename' (easier on the mechanism loads)
*/

// main experimental parameters and parameter initialization

max_branch_levels = 7 //  maximum nr. of branch levels (from 0 to 6)

branch_levels = 1 // nr of branch levels in dendrite (if 0, only apical dendrite is present)

begin_branch_diam = 2.5 // diameter micrometer of apical dendrite
exponent = 3/2 // diameter rule exponent
max_seg_length = 30 // max length of segments (micrometer)


//TotalPump (total density of pump sites on the cell membrane (mol/cm2))
TotalPump_cadp = 0.3e-12

// TotalBuffer (total amount of buffer present)
TotalBuffer_cadp = 0.000

v_init = -70 // initial membrane potential mV (same as rest potential

NDEND = 0 // number of dendrites
NDEND2 = 0 // will be 1 if branch_levels is set to zero 
// in that case an unconnected piece of dendrite is added to the model (treat exceptional case)


// compute the number of dendritic segments (apart from the apical segment)
for i=1, branch_levels {
    NDEND=NDEND+(2^i)
}

NDEND2 = NDEND

// if branch_levels == 0, (NDEND==0) then treat exceptional case (add an unconnected piece of dendrite)
if (NDEND<1){
    NDEND2=1 // one disconnected segment will be created
}

objref leveldiam
leveldiam = new Vector(branch_levels+1,0)

objref levellength
levellength = new Matrix(max_branch_levels+1,max_branch_levels+1)

// Topology
create soma, apical, dend[NDEND2], hillock
access soma

connect soma(0), hillock(1)
connect apical(0), soma(1)

// connect dendritic segments if branch_levels
if (branch_levels>0) {
    for i = 0, branch_levels-1{
        start_dend = 0

        for k = 0, i{
            start_dend=start_dend+2^k
        }

        start_dend=start_dend-1

        for p=0, (2^(i+1))-1{
            if (i<1){
                connect dend[start_dend+p](0),apical(1)

            } else{
                connect dend[start_dend+p](0),dend[(start_dend-(2^i))+int(p/2)](1)
            }
        }
    }
}


// Geometry

path_length_dendr = 250 // path length along dendrite

sc = 0.5 // scaling factor of length between branch levels

// data array of dendritic compartment lengths (branch level dependent)
levellength.x[0][0] = path_length_dendr
levellength.x[1][0] = path_length_dendr - (path_length_dendr*sc)
levellength.x[1][1] = path_length_dendr*sc
levellength.x[2][0] = path_length_dendr - (path_length_dendr*sc)
levellength.x[2][1] = path_length_dendr*sc - (path_length_dendr*(sc^2))
levellength.x[2][2] = path_length_dendr*(sc^2)
levellength.x[3][0] = path_length_dendr - (path_length_dendr*sc)
levellength.x[3][1] = path_length_dendr*sc - (path_length_dendr*(sc^2))
levellength.x[3][2] = (path_length_dendr*(sc^2))-(path_length_dendr*(sc^3))
levellength.x[3][3] = path_length_dendr*(sc^3)
levellength.x[4][0] = path_length_dendr - (path_length_dendr*sc)
levellength.x[4][1] = path_length_dendr*sc - (path_length_dendr*(sc^2))
levellength.x[4][2] = (path_length_dendr*(sc^2)) - (path_length_dendr*(sc^3))
levellength.x[4][3] = (path_length_dendr*(sc^3))- (path_length_dendr*(sc^4))
levellength.x[4][4] = (path_length_dendr*(sc^4))
levellength.x[5][0] = path_length_dendr - (path_length_dendr*sc)
levellength.x[5][1] = path_length_dendr*sc - (path_length_dendr*(sc^2))
levellength.x[5][2] = (path_length_dendr*(sc^2)) - (path_length_dendr*(sc^3))
levellength.x[5][3] = (path_length_dendr*(sc^3))- (path_length_dendr*(sc^4))
levellength.x[5][4] = (path_length_dendr*(sc^4))- (path_length_dendr*(sc^5))
levellength.x[5][5] = (path_length_dendr*(sc^5))
levellength.x[6][0] = path_length_dendr - (path_length_dendr*sc)
levellength.x[6][1] = path_length_dendr*sc - (path_length_dendr*(sc^2))
levellength.x[6][2] = (path_length_dendr*(sc^2)) - (path_length_dendr*(sc^3))
levellength.x[6][3] = (path_length_dendr*(sc^3))- (path_length_dendr*(sc^4))
levellength.x[6][4] = (path_length_dendr*(sc^4))- (path_length_dendr*(sc^5))
levellength.x[6][5] = (path_length_dendr*(sc^5)) - (path_length_dendr*(sc^6))
levellength.x[6][6] = (path_length_dendr*(sc^6))
levellength.x[7][0] = path_length_dendr - (path_length_dendr*sc)
levellength.x[7][1] = path_length_dendr*sc - (path_length_dendr*(sc^2))
levellength.x[7][2] = (path_length_dendr*(sc^2)) - (path_length_dendr*(sc^3))
levellength.x[7][3] = (path_length_dendr*(sc^3))- (path_length_dendr*(sc^4))
levellength.x[7][4] = (path_length_dendr*(sc^4))- (path_length_dendr*(sc^5))
levellength.x[7][5] = (path_length_dendr*(sc^5)) - (path_length_dendr*(sc^6))
levellength.x[7][6] = (path_length_dendr*(sc^6)) - (path_length_dendr*(sc^7))
levellength.x[7][7] = (path_length_dendr*(sc^7))

nr_end_branches = 2^(branch_levels) // number of end points / branches

// compute total length of dendritic tree
total_length = 0 // initial value for sum of total_length
for i = 0, branch_levels{
	total_length = total_length+(levellength.x[branch_levels][i]*(2^i))
}
	
//compute DCI of dendritic trees
DCI = (nr_end_branches+(nr_end_branches*branch_levels))*total_length

print "dci"
print DCI
print "total length"
print total_length

// the Rall value to keep constant
diameter_constant = begin_branch_diam^exponent

// compute diameters for various branch levels (apply Rall rule)
for i=0, branch_levels {
    leveldiam.x[i] = (diameter_constant/(2^i))^(1/exponent)
    print leveldiam.x[i]
}

dend_L = path_length_dendr/(branch_levels+1)

//apical geometry
apical{
    L = levellength.x[branch_levels][0]
    diam = leveldiam.x[0]
    nseg = int(levellength.x[branch_levels][0]/max_seg_length)+1
}

//set dendritic geometry (diameter, length, nr of segments
if (branch_levels>0){

    for i=1, branch_levels{
        start_dend = 0


        for k = 0, i-1 {
        start_dend=start_dend+2^k
        }

        start_dend=start_dend-1

        for p=0, (2^(i))-1{

            dend[start_dend+p]{

                diam = leveldiam.x[i]
                L = levellength.x[branch_levels][i]
                nseg = int(levellength.x[branch_levels][i]/max_seg_length)+1

            }
        }
    }
}

soma_L = 10 // somatic length (micrometer)
soma_diam = 15 // soma diameter (micrometer)

hillock_L = 2 // axon hillock length (micrometer)
hillock_diam = 2 // axon hillock diameter (micrometer)

// set soma geometric properties
soma{
    L = soma_L
    diam = soma_diam
    nseg = 1
}

// set axon hillock geometric properties
hillock{
    L = hillock_L
    diam = hillock_diam
}



// Biophysics
ra        = 150 // axial resistance
global_ra = ra
rm        = 30000  // 30 kilo ohm
c_m       = 0.75 // micro farad

// passive apical properties
apical {
    insert pas
    Ra = ra
    cm = c_m
    g_pas = 1/rm
    e_pas = v_init
}

// passive soma properties
soma {
    insert pas
    Ra = ra
    cm = c_m
    g_pas = 1/rm
    e_pas = v_init
}


celsius = 32 // degrees celsius of experiment

Ek = -90 // reversal potential potassium
Ena = 50 // reversal potential sodium


// dendritic conductances 
gna = 1 
gca = 2  
gkca = 10
gkm = 0.1 

// soma conductances
gkv_soma= 800
gna_soma= 300
gca_soma = 0.1
gkca_soma = 0
gkm_soma = 0

//hillock conductances
gkv_hillock = 2000
gna_hillock = 30000


// somatic active channels

soma{
    insert na gbar_na = gna_soma
    insert kv gbar_kv = gkv_soma
    insert ca gbar_ca = gca_soma
    insert cadp
}


// apical active channels
apical{
    insert na	 gbar_na = gna
    insert km    gbar_km  = gkm
    insert ca    gbar_ca = gca
    insert cadp
	insert kca   gbar_kca = gkca

}
  
// axon hilock active channels  
hillock{
	insert na gbar_na = gna_hillock
	insert kv gbar_kv = gkv_hillock
}

// insert passive and active properties in dendritic compartments
for i=0, NDEND-1{
    dend[i]{     
	    insert pas
        Ra = ra
        cm = c_m
        g_pas = 1/rm
        e_pas = v_init
	    insert na	 gbar_na = gna
        insert km    gbar_km  = gkm
        insert ca    gbar_ca = gca
        insert cadp
	    insert kca   gbar_kca = gkca
    }
}

// set general extracellular properties for calcium
forall if(ismembrane("ca_ion")) {
    eca = 140
    ion_style("ca_ion",0,1,0,0,0)
    vshift_ca = 0
}



// Instrumentation

// current input
objref stim
soma stim = new IClamp(0.5)
stim.amp = .02
stim.del = 50
stim.dur = 1000

//Simulation control
tstop = 2000
steps_per_ms = 40
dt = 0.025
nr_cells=int(tstop/dt)+2

// adjust output data size to number of branch levels
if(branch_levels<5){
    nvars = 5 // no thin dendrite calcium dynamics
}else{
    nvars = 6
}
double data[nvars][nr_cells] // create empty array for data storage


// Plotting

// graph of somatic membrane potential
objref g
g = new Graph()
g.size(0,1050,-70,60)
g.addvar("soma.v(0.5)",1,1,0.6,0.9,2)
graphList[1].append(g)
g.save_name("somatic voltage against time")

// Experimental control
proc initialize() {
    t = 0
    cnt = 0
    finitialize(v_init)
    fcurrent()
}

proc integrate() {
    g.begin()
    while (t<tstop) {
        // next section write to the 'data' array
        if(branch_levels<5){
            data[0][cnt]= t
            data[1][cnt] = soma.v(.5)
            data[2][cnt] = soma.cai(.5)
            data[3][cnt] = apical.cai(0.5)
            data[4][cnt] = dend[0].cai(.5) // 1st branch level calcium readout
        }else{
            data[0][cnt]= t
            data[1][cnt] = soma.v(.5)
            data[2][cnt] = soma.cai(.5)
            data[3][cnt] = apical.cai(0.5)
            data[4][cnt] = dend[0].cai(.5)
            data[5][cnt] = dend[14].cai(.5) // 4th branch level calcium readout
        }
        cnt=cnt+1
        g.plot(t)
        fadvance()

    }
    g.flush()
}

// function runs experiment
proc go() {
    initialize()
    integrate()
    finalize()
}

// finalize writes data to file
proc finalize(){

    strdef filename, prefix, postfix
    prefix = "altered_complexity_branch_levels_ratio"
    postfix="txt"
    num = branch_levels
    num2 = exponent
    sprint(filename, "%s_%d_%f.%s", prefix, num, num2,postfix)
    wopen(filename)

    // write data to text file
    for i=0,nr_cells-1{
        if(branch_levels<5){
	        fprint("%f\t %f\t %f\t %f\t %f\n",data[0][i],data[1][i],data[2][i],data[3][i],data[4][i])
        }else{
	        fprint("%f\t %f\t %f\t %f\t %f\t %f\n",data[0][i],data[1][i],data[2][i],data[3][i],data[4][i],data[5][i])
        }
    }
    wopen()
}

go() // run experiment