/* fig6a3-4.hoc Driver Program for Midbrain Dopaminergic Neuron model */

/* The following command loads the standard run library functions into
 * the simulator.  These library functions are used to initialize the
 * simulator, begin and end simulations, and plot variables during a
 * simulation. */
timecon = 0
back= 1                         /* 1 switch for running in the background */
sakmann=0
if(!back)load_file("nrngui.hoc")
verbose=0
v_init = -60
vcseries = 0
clamp = 0                     /* switch for voltage clamp*/
restart = 0                      /* switch for initializing from state.old */
tstart = 0
if(vcseries) tstop = 800 else tstop = 10000
if(clamp && !vcseries) tstop= 600
/*if(timecon) tstop = 100*/
if(timecon) tstop = 800
nsyn = 1			/* The number of synapses */
soma_len= 25.0 /*25*/
soma_diam= 15.0
ndend=4
nbranch=2
dend_diam= 1.5  /* 1.5 need to make dendrite taper */
prox_diam= 3.0  /* 3.0 need to make dendrite taper */
taper = 1.0
dend_len = 500.0
prox_len = 150.0
Dtmax =   1.0  
Dt =   1.00
if(timecon) Dtmax = 1.0 
dt = 5e-4  /*0.1 5e-9 */
if(timecon) dt = 0.01
nainit=  4.000
vsolder=v_init
vdolder=v_init
vsold=v_init
vdold=v_init

create dummy,soma, prox[ndend], dend[nbranch*ndend]	/* Create all of the cell's sections */
global_ra = 400
forall Ra = global_ra
global_cm = 1.0
forall cm = global_cm
g_celsius = 35
celsius = g_celsius
forall ion_style("na_ion", 2,2,0,0,0)
access soma			/* All default references are to soma */

objectvar syn[nsyn]		/* Declare the object variables for
				 * the three synapses */
objectvar stim
objectvar vc
objref cvode
cvode = new CVode(1) /* 0 for clamp*/
x= cvode.active(1)

proc init_cell() {
/* First set all of the dimensions and insert the channels into each
 * section. */
    dummy{
      L=1
      d=1
    }
    soma {
	nseg = 1		/* # of compartments */
	diam = soma_diam		/* Set the diameter of the soma (in um) */
	L = soma_len			/* Set the length of the soma (in um) */
        
	{insert nadifl D_nadifl = 0.60 nainit_nadifl = 2.4 f_nadifl = 4.0}
        {insert hh3 gnabar_hh3 = 5500.0e-6 gkabar_hh3 = 100.0e-6 gkhhbar_hh3 = 1000.0e-6
          miv_hh3 = 44.6 htv1_hh3 = 39.0 hiv_hh3 = 66.8 htv2_hh3 = 59.0}
        {insert pump ipumpmax_pump = 0.0036 km_pump = 10.0  }
        {insert leak gkbar_leak = 18.0e-6  gnabar_leak = 9.5e-6  
          ggabaa_leak = 1450.0e-6 gcabar_leak = 0.6e-6 }
        {insert cadifus fCa_cadifus = 0.005 cainit_cadifus = 0.000116 }
       {insert calcha gcatbar_calcha = 1044.0e-6 gcanbar_calcha = 171.0e-6 
          gcalbar_calcha = 216.0e-6 gcahvabar_calcha = 0.0e-6}
        {insert gkca gkbar_gkca = 900.0e-6 /*0.0e-6*/
	 ek_gkca = -100 km_gkca = 0.00019
          n_gkca = 4.0}
        {insert capump icapumpmax_capump = 0.0312 km_capump = 0.0005}
        if (!clamp) {
        istim = 0.0 /* 0.07 nA=70pA*/
        stim = new MyIClamp(.5)
        {stim.del = 0 stim.dur = 200000 stim.amp = 0.0
         stim.amp2 = 0}
if(timecon) stim.dur = 600  /*(ms)*/
if (sakmann) stim.amp2 = -0.050
if(timecon) stim.amp2 = -0.020 
if(timecon) stim.del = 100
     } else {
	vc = new SEClamp(.5)
        if(!vcseries) vc.amp1 = -60.0 else vc.amp1 = -130
	vc.dur1 =  50
	vc.amp2 = -30.0
	vc.dur2 = 350
	vc.amp3 = -60.0
	vc.dur3 = 201
    }

	pt3dclear()
	pt3dadd(0,0,0,soma_diam)
	pt3dadd(soma_len,0,0,soma_diam)

    }
 
    for i = 0, ndend-1 prox[i] {
	nseg = 1	/* # of compartments */
	L = prox_len			/* Set the length of the dendrite (in um) */
     
	{insert nadifl D_nadifl = 0.60  nainit_nadifl = 3.0  f_nadifl = 1.0}
	{insert hh3 gnabar_hh3 = 5500.0e-6 gkabar_hh3 = 300.0e-6 gkhhbar_hh3 = 1000.0e-6
          miv_hh3 = 34.6 htv1_hh3 = 29.0 hiv_hh3 = 56.8 htv2_hh3 = 49.0}
	{insert pump ipumpmax_pump = 0.0072  km_pump = 10.0}
	{insert leak gkbar_leak = 18.0e-6 gnabar_leak = 9.5e-6 
          ggabaa_leak = 145.0e-6 gcabar_leak = 0.0e-6 }
        {insert nmda Pbar_nmda = 2.5e-6 }
        {insert ampa gampa_ampa = 0.0e-6 gampak_ampa = 0.0e-6 ek_ampa = -100.0 ratio = 1}
        
forall cm = global_cm
forall Ra = global_ra
g_celsius = 35
    }
 
    for i = 0, nbranch*ndend - 1 dend[i] {
	nseg = 1	/* # of compartments */
	L = dend_len - prox_len 			/* Set the length of the dendrite (in um) */
        
	{insert nadifl D_nadifl = 0.60 nainit_nadifl = 4.59 f_nadifl = 1.0}
	{insert hh3 gnabar_hh3 = 5500.0e-6 gkabar_hh3 = 1000.0e-6 gkhhbar_hh3 = 1000e-6
           miv_hh3 = 26.6 htv1_hh3 = 21 hiv_hh3 = 48.8 htv2_hh3 = 41.0}
        {insert pump ipumpmax_pump = 0.009  km_pump = 10.0}
        {insert leak gkbar_leak = 18.0e-6 gnabar_leak = 9.5e-6  
           ggabaa_leak = 145.0e-6 gcabar_leak = 0.0e-6 }
        {insert nmda Pbar_nmda = 2.5e-6}
	{insert ampa  gampa_ampa = 0.0e-6 gampak_ampa = 0.0e-6 ek_ampa = -100.0 ratio = 1}
    

forall cm = global_cm
forall Ra = global_ra
g_celsius = 35
    }

/* Next we construct the topology by connecting each of the sections
 * together. */
    connect soma(0),dummy(1)
    connect prox[0](0),soma(0)
    connect prox[1](0), soma(1)
    connect prox[2](0), soma(0.5)
    connect prox[3](0), soma(0.5) 
    for j = 0, ndend-1 {
    for i = 0, nbranch-1 {
    connect dend[j+i*ndend](0), prox[j](1)}}

        dummy{
	pt3dclear()
	pt3dadd(0.5*soma_len,0,1,1)
	pt3dadd(0.5*soma_len,0,0,1)
               }
        prox[0]{
	pt3dclear() 
	pt3dadd(0,0,0,prox_diam)
	pt3dadd(-prox_len,0,0,prox_diam)
               }
        dend[0]{
	pt3dclear() 
	pt3dadd(-prox_len,0,0,dend_diam)
	pt3dadd(-dend_len,0,0,dend_diam/taper)
               }
        if(nbranch==2) dend[4]{
	pt3dclear() 
	pt3dadd(-prox_len,0,0,dend_diam)
	pt3dadd(-prox_len,dend_len-prox_len,0,dend_diam/taper)
               }
        prox[1]{
	pt3dclear() 
	pt3dadd(soma_len,0,0,prox_diam)
	pt3dadd(soma_len+prox_len,0,0,prox_diam)
               }
        dend[1]{
	pt3dclear() 
	pt3dadd(soma_len+prox_len,0,0,dend_diam)
	pt3dadd(soma_len+dend_len,0,0,dend_diam/taper)
               }
        if(nbranch==2) dend[5]{
	pt3dclear() 
	pt3dadd(soma_len+prox_len,0,0,dend_diam)
	pt3dadd(soma_len+prox_len,prox_len-dend_len,0,dend_diam/taper)
               }
        prox[2]{
	pt3dclear() 
	pt3dadd(0,soma_diam/2,0,prox_diam)
	pt3dadd(0,(soma_diam/2 + prox_len),0,prox_diam)
               }
        dend[2]{
	pt3dclear() 
	pt3dadd(0,(soma_diam/2 + prox_len),0,dend_diam)
	pt3dadd(0,soma_diam/2 + dend_len,0,dend_diam/taper)
               }
        if(nbranch==2) dend[6]{
	pt3dclear() 
	pt3dadd(0,(soma_diam/2 + prox_len),0,dend_diam)
	pt3dadd(dend_len-prox_len,soma_diam/2 + prox_len,0,dend_diam/taper)
               }
        prox[3]{
	pt3dclear() 
	pt3dadd(0,-soma_diam/2,0,prox_diam)
	pt3dadd(0,(-soma_diam/2 - prox_len),0,prox_diam)
               }
        dend[3]{
	pt3dclear() 
	pt3dadd(0,-soma_diam/2 - prox_len,0,dend_diam)
	pt3dadd(0,-soma_diam/2 - dend_len,0,dend_diam/taper)
               }
        if(nbranch==2) dend[7]{
	pt3dclear() 
	pt3dadd(0,-soma_diam/2 - prox_len,0,dend_diam)
	pt3dadd(prox_len-dend_len,-soma_diam/2 - prox_len,0,dend_diam/taper)
               }
}


topology()

init_cell()			/* Call the function to initialize our
				 * cell. */
objref ss,f1,f2
ss = new SaveState()
f1 = new File()
f2 = new File()



proc init() {local i
if(!restart){
             finitialize(v_init)
             fcurrent()
             }
if(restart){ 
           f1.ropen("state.old")
           ss.fread(f1)
           f1.close
           finitialize(v_init)
           ss.restore()
           t=tstart
           if(cvode.active()){
           cvode.re_init()

} else {
           fcurrent()
        }
 frecord_init()
           }
        t = tstart
}




/*proc continuerun() {local rt
	eventcount=0
	eventslow=1
	stoprun = 0
        startsw()
	while(t < $1 && stoprun == 0) {
		step()
		rt = stopsw()
		if (rt > realtime) {
			realtime = rt
			fastflushPlot()
			doNotify()
			if (realtime == 2 && eventcount > 50) {
				eventslow = int(eventcount/50) + 1
			}
			eventcount = 0
		}else{
			eventcount = eventcount + 1
			if ((eventcount%eventslow) == 0) {
				doEvents()
			}
		}
	}
f2.wopen("state.new")
ss.save()
ss.fwrite(f2)
f2.close
flushPlot()
}*/
init()
/*startsw()
flushPlot()*/
if(back){
if(!clamp ||verbose) {print t,soma.v(0.5),dend[1].v(0.99),soma.nai(0.5), dend[1].nai(0.99),soma.cai(0.5),dend[1].inapump_pump(0.99), dend[1].inmda_nmda(0.99)}
if(clamp && !vcseries && !verbose) print t,vc.i,soma.cai(0.5),soma.v(0.5)
if(vcseries) j= 10 else j = 0
for i = 0, j { if (vcseries) vc.amp1 = vc.amp1 + 10}
if(vcseries) init()
next = t + Dt
flag1=0
flag2=0
hold = 0
while (t<tstop){
          vsolder=vsold
          vdolder=vdold
          vsold=soma.v(0.5)
          vdold=dend[1].v(0.99)
          fadvance()
           if(!clamp||verbose){
          if((vsolder<vsold &&soma.v(0.5) <vsold)||(vsolder>vsold && soma.v(0.5)>vsold)||(vdolder<vdold &&dend[1].v(0.99) <vdold)|| (vdolder>vdold && dend[1].v(0.99)>vdold)) {
          vsolder=soma.v(0.5)
          vdolder=dend[1].v(0.99)
                           print t,soma.v(0.5),dend[1].v(0.99),soma.nai(0.5), dend[1].nai(0.99),soma.cai(0.5), dend[1].inapump_pump(0.99),dend[1].inmda_nmda(0.99)
               next = t + Dt
                           flag2=1
                           hold=dt
          soma.v(0.5)=vsolder
          dend[1].v(0.99)=vdolder
                            }
                            }   
               if(t>=next){
               Dt = 100*dt
               if(Dt>Dtmax) Dt = Dtmax
               if (Dt<.1) Dt = .1
               next = t + Dt
         if(!clamp||verbose) print t,soma.v(0.5),dend[1].v(0.99),soma.nai(0.5), dend[1].nai(0.99),soma.cai(0.5), dend[1].inapump_pump(0.99),dend[1].inmda_nmda(0.99)
if(clamp && !vcseries && !verbose) print t,vc.i,soma.cai(0.5),soma.v(0.5)
                          }
                 }
if(vcseries) print vc.amp1,vc.i
f2.wopen("state.new")
ss.save()
ss.fwrite(f2)
f2.close
} else{ if(clamp) {xopen("clamp.ses")} else {xopen("new.ses")}
                   forall Ra = global_ra
                   forall cm = global_cm
                   prox{ forall cm = global_cm}
                   dend{ forall cm = global_cm}
                   celsius = g_celsius
       }