// this hoc file represents a convertion to NEURON from the Matlab model originally published by Renaud Jolivet (2015)
// this hoc file drives coggan_keller_joliver_orig_Modeldb.model

load_file("nrngui.hoc")
create soma
access soma

insert glia_2013
time =0
tstop=7500000.0//ms
tend=tstop
tfinal=tstop

tstim=5000000.0 //ms
tstimend=00006000000//ms
tBloodEnd=60000 //ms

dt=.1//ms

soma.Isyn_glia_2013=0

soma.F0_glia_2013=0.012

finitialize()

gRate=0.0002
F0=0.012
lag=1000
t1=0

//add the glutamatergic stimulation
objref output_file
output_file = new File()
output_file.wopen("output_file.dat")

gbar    = 7.8e-06
Ne=240

skipstep=0
oldy=0
for i=0,tstop/( dt) {
  fadvance()
  if (time>tstim){
    if (time<tstimend){

      soma.Nag_glia_2013=soma.Nag_glia_2013+gRate*dt

    finf    = 2.5  // Hz
    f0      = 23   // Hz
    tau     = 2500 //msec

   //soma.Isyn_glia_2013= Ne*gbar*(finf+(f0-finf)*exp(-time/tau))
soma.Isyn_glia_2013= 0

    }
  }
  if (time>tstimend){
    soma.Isyn_glia_2013=0
  }
  if (skipstep>100/dt){


//////////////////

	
//1
output_file.printf("%.10g ",time/1000)
//2
output_file.printf("%.10g ",soma.Nan_glia_2013(0.5))
//3
output_file.printf("%.10g ",soma.Nag_glia_2013(0.5))
//4
output_file.printf("%.10g ",soma.GLCn_glia_2013(0.5))
//5
output_file.printf("%.10g ",soma.GLCg_glia_2013(0.5))
//6
output_file.printf("%.10g ",soma.GAPn_glia_2013(0.5))
//7
output_file.printf("%.10g ",soma.GAPg_glia_2013(0.5))
//8
output_file.printf("%.10g ",soma.PEPn_glia_2013(0.5))
//9
output_file.printf("%.10g ",soma.PEPg_glia_2013(0.5))
//10
output_file.printf("%.10g ",soma.PYRn_glia_2013(0.5))
//11
output_file.printf("%.10g ",soma.PYRg_glia_2013(0.5))
//12
output_file.printf("%.10g ",soma.LACn_glia_2013(0.5))
//13
output_file.printf("%.10g ",soma.LACg_glia_2013(0.5))
//14
output_file.printf("%.10g ",soma.NADHcyton_glia_2013(0.5))
//15
output_file.printf("%.10g ",soma.NADHcytog_glia_2013(0.5))
//16
output_file.printf("%.10g ",soma.ATPn_glia_2013(0.5))
//17
output_file.printf("%.10g ",soma.ATPg_glia_2013(0.5))
//18
output_file.printf("%.10g ",soma.PCrn_glia_2013(0.5))
//19
output_file.printf("%.10g ",soma.PCrg_glia_2013(0.5))
//20
output_file.printf("%.10g ",soma.O2n_glia_2013(0.5))
//21
output_file.printf("%.10g ",soma.O2g_glia_2013(0.5))
//22
output_file.printf("%.10g ",soma.O2c_glia_2013(0.5))
//23
output_file.printf("%.10g ",soma.GLCc_glia_2013(0.5))
//24
output_file.printf("%.10g ",soma.LACc_glia_2013(0.5))
//25
output_file.printf("%.10g ",soma.Vv_glia_2013(0.5))
//26
output_file.printf("%.10g ",soma.Hb_glia_2013(0.5))
//27
output_file.printf("%.10g ",soma.GLCe_glia_2013(0.5))
//28
output_file.printf("%.10g ",soma.LACe_glia_2013(0.5))
//29
output_file.printf("%.10g ",soma.Ca_glia_2013(0.5))
//30
output_file.printf("%.10g ",soma.NADHmiton_glia_2013(0.5))
//31
output_file.printf("%.10g ",soma.NADHmitog_glia_2013(0.5))
//32
output_file.printf("%.10g ",soma.Rminusn_glia_2013(0.5))
//33
output_file.printf("%.10g ",soma.Rplusn_glia_2013(0.5))
//34
output_file.printf("%.10g ",soma.Rminusg_glia_2013(0.5))
//35
output_file.printf("%.10g ",soma.Rplusg_glia_2013(0.5))
//36
output_file.printf("%.10g ",soma.JleakNag_glia_2013(0.5))
//37
output_file.printf("%.10g ",soma.JleakNan_glia_2013(0.5))
//38
output_file.printf("%.10g ",soma.Jpumpg_glia_2013(0.5))
//39
output_file.printf("%.10g ",soma.Jpumpn_glia_2013(0.5))
//40
output_file.printf("%.10g ",soma.JGLCce_glia_2013(0.5))
//41
output_file.printf("%.10g ",soma.JGLCcg_glia_2013(0.5))
//42
output_file.printf("%.10g ",soma.JGLCeg_glia_2013(0.5))
//43
output_file.printf("%.10g ",soma.JGLCen_glia_2013(0.5))
//44
output_file.printf("%.10g ",soma.JHKPFKn_glia_2013(0.5))
//45
output_file.printf("%.10g ",soma.JHKPFKg_glia_2013(0.5))
//46
output_file.printf("%.10g ",soma.ADPn_glia_2013(0.5))
//47
output_file.printf("%.10g ",soma.ADPg_glia_2013(0.5))
//48
output_file.printf("%.10g ",soma.JPGKn_glia_2013(0.5))
//49
output_file.printf("%.10g ",soma.JPGKg_glia_2013(0.5))
//50
output_file.printf("%.10g ",soma.JPKn_glia_2013(0.5))
//51
output_file.printf("%.10g ",soma.JPKg_glia_2013(0.5))
//52
output_file.printf("%.10g ",soma.JLDHn_glia_2013(0.5))
//53
output_file.printf("%.10g ",soma.JLDHg_glia_2013(0.5))
//54
output_file.printf("%.10g ",soma.JLACec_glia_2013(0.5))
//55
output_file.printf("%.10g ",soma.JLACcg_glia_2013(0.5))
//56
output_file.printf("%.10g ",soma.JLACeg_glia_2013(0.5))
//57
output_file.printf("%.10g ",soma.JLACen_glia_2013(0.5))
//58
output_file.printf("%.10g ",soma.Jmitoing_glia_2013(0.5))
//59
output_file.printf("%.10g ",soma.Jmitoinn_glia_2013(0.5))
//60
output_file.printf("%.10g ",soma.Jmitooutn_glia_2013(0.5))
//61
output_file.printf("%.10g ",soma.Jmitooutg_glia_2013(0.5))
//62
output_file.printf("%.10g ",soma.Jshuttlen_glia_2013(0.5))

// additions below

//63/ 109
output_file.printf("%.10g ",soma.psin_glia_2013(0.5))
//64/ 110
output_file.printf("%.10g ",soma.Jstimn_glia_2013(0.5))
//65 /111
output_file.printf("%.10g ",soma.Jstimg_glia_2013(0.5))
//66 /112
output_file.printf("%.10g ",soma.ICa_glia_2013(0.5))
//67 /113
output_file.printf("%.10g ",soma.INa_glia_2013(0.5))
//68 /114
output_file.printf("%.10g ",soma.IK_glia_2013(0.5))



//69 /120
output_file.printf("%.10g ",soma.JLDHg_glia_2013(0.5))
//70 /121
output_file.printf("%.10g ",soma.JLACge_glia_2013(0.5))
//71/ 122
output_file.printf("%.10g ",soma.JLACgc_glia_2013(0.5))
//72 /123
output_file.printf("%.10g ",soma.NADHcytog_glia_2013(0.5))
//73 /124
output_file.printf("%.10g ",soma.LACe_glia_2013(0.5))
//74 /125
output_file.printf("%.10g ",soma.TmaxLACeg_glia_2013(0.5))
//75 /126
output_file.printf("%.10g ",soma.TmaxLACge_glia_2013(0.5))
//76
output_file.printf("%.10g ",soma.JCKg_glia_2013(0.5))
//77
output_file.printf("%.10g ",soma.JCKn_glia_2013(0.5))

//78
output_file.printf("%.10g ",soma.dAMPdATPn_glia_2013(0.5))

//79
output_file.printf("%.10g ",soma.dAMPdATPg_glia_2013(0.5))
//80
output_file.printf("%.10g ", soma.un_glia_2013(0.5))
//81  
output_file.printf("%.10g ", soma.ug_glia_2013(0.5))
//82 
output_file.printf("%.10g ", soma.Isyn_glia_2013(0.5))
//83
//output_file.printf("%.10g ", soma.vnstim_glia_2013(0.5))
//84
//output_file.printf("%.10g ", soma.vgstim_glia_2013(0.5))

//83
output_file.printf("%.10g ",soma.Fin_glia_2013(0.5))
//84
output_file.printf("%.10g \n",soma.Finprime_glia_2013(0.5))


/////////////

    skipstep=0

  }
  skipstep=skipstep+1

 if (1>0) {
  //if (1<0){
    //put in cortical blood flow 

    lastF0 = F0*( 1.1+1.5*(exp(-(tBloodEnd-lag)/5000)-exp(-(tBloodEnd-lag)/2000)) ) //builtin1(tBloodEnd-lag);
    taurise= 1000 * 0.1*1/( 1.5*(-1/5*exp(-0/5)+1/2*exp(-0/2))) //builtin2(0);

    if ((time <= (tBloodEnd+t1)) && (time >= (t1+lag))){
            y   =       F0* (1.1+1.5*(exp(-(time-lag)/5000)-exp(-(time-lag)/2000))) //builtin1(time-lag);
    }else{
      if (time > (tBloodEnd)){
          y   =   F0+(lastF0-F0)*exp(-(time-tBloodEnd)/5000)
      } else{
           y   =   F0*(1+0.1*exp((time-lag)/taurise)) //builtin3(t,lag,taurise))
      }
  }
  soma.Fin_glia_2013=y
  if (i>0){
    soma.Finprime_glia_2013=(y-oldy)/dt
  
  }
  oldy=y
}
//take these next two lines out when running real sim with blood flow modulation
soma.Finprime_glia_2013=0
soma.Fin_glia_2013=0.012

  time=time+ dt
}


output_file.close()

output_file.wopen("output_file2.dat")
 

output_file.printf("Nan= %.6g\n", soma.Nan_glia_2013(0.5))
output_file.printf("Nag= %.6g\n", soma.Nag_glia_2013(0.5))
output_file.printf("GLCn= %.6g\n", soma.GLCn_glia_2013(0.5))
output_file.printf("GLCg= %.6g\n", soma.GLCg_glia_2013(0.5))
output_file.printf("GAPn= %.6g\n", soma.GAPn_glia_2013(0.5))
output_file.printf("GAPg= %.6g\n", soma.GAPg_glia_2013(0.5))
output_file.printf("PEPn= %.6g\n", soma.PEPn_glia_2013(0.5))
output_file.printf("PEPg= %.6g\n", soma.PEPg_glia_2013(0.5))
output_file.printf("PYRn= %.6g\n", soma.PYRn_glia_2013(0.5))
output_file.printf("PYRg= %.6g\n", soma.PYRg_glia_2013(0.5))
output_file.printf("LACn= %.6g\n", soma.LACn_glia_2013(0.5))
output_file.printf("LACg= %.6g\n", soma.LACg_glia_2013(0.5))
output_file.printf("NADHcyton= %.6g\n", soma.NADHcyton_glia_2013(0.5))
output_file.printf("NADHcytog= %.6g\n", soma.NADHcytog_glia_2013(0.5))
output_file.printf("ATPn= %.6g\n", soma.ATPn_glia_2013(0.5))
output_file.printf("ATPg= %.6g\n", soma.ATPg_glia_2013(0.5))
output_file.printf("PCrn= %.6g\n", soma.PCrn_glia_2013(0.5))
output_file.printf("PCrg= %.6g\n", soma.PCrg_glia_2013(0.5))
output_file.printf("O2n= %.6g\n", soma.O2n_glia_2013(0.5))
output_file.printf("O2g= %.6g\n", soma.O2g_glia_2013(0.5))
output_file.printf("O2c= %.6g\n", soma.O2c_glia_2013(0.5))
output_file.printf("GLCc= %.6g\n", soma.GLCc_glia_2013(0.5))
output_file.printf("LACc= %.6g\n", soma.LACc_glia_2013(0.5))
output_file.printf("Vv= %.6g\n", soma.Vv_glia_2013(0.5))
output_file.printf("Vvp= %.6g\n", soma.Vvp_glia_2013(0.5))
output_file.printf("Hb= %.6g\n", soma.Hb_glia_2013(0.5))
output_file.printf("GLCe= %.6g\n", soma.GLCe_glia_2013(0.5))
output_file.printf("LACe= %.6g\n", soma.LACe_glia_2013(0.5))
output_file.printf("psin= %.6g\n", soma.psin_glia_2013(0.5))
output_file.printf("h= %.6g\n", soma.h_glia_2013(0.5))
output_file.printf("n= %.6g\n", soma.n_glia_2013(0.5))
output_file.printf("Ca= %.6g\n", soma.Ca_glia_2013(0.5))
output_file.printf("NADHmiton= %.6g\n", soma.NADHmiton_glia_2013(0.5))
output_file.printf("NADHmitog= %.6g\n", soma.NADHmitog_glia_2013(0.5))
output_file.printf("Fout= %.6g\n", soma.Fout_glia_2013(0.5))
output_file.printf("\n")


output_file.printf("ADPn= %.6g\n", soma.ADPn_glia_2013(0.5))
output_file.printf("ADPg= %.6g\n", soma.ADPg_glia_2013(0.5))
  output_file.printf("Rminusn = %.6g\n", soma.Rminusn_glia_2013(0.5))

  output_file.printf("Rplusn = %.6g\n", soma.Rplusn_glia_2013(0.5))

  output_file.printf("Rminusg = %.6g\n", soma.Rminusg_glia_2013(0.5))

  output_file.printf("Rplusg = %.6g\n", soma.Rplusg_glia_2013(0.5))

  output_file.printf("JleakNag = %.6g\n", soma.JleakNag_glia_2013(0.5))

  output_file.printf("JleakNan = %.6g\n", soma.JleakNan_glia_2013(0.5))

  output_file.printf("Jpumpg = %.6g\n", soma.Jpumpg_glia_2013(0.5))

  output_file.printf("Jpumpn = %.6g\n", soma.Jpumpn_glia_2013(0.5))

  output_file.printf("JGLCce = %.6g\n", soma.JGLCce_glia_2013(0.5))

  output_file.printf("JGLCcg = %.6g\n", soma.JGLCcg_glia_2013(0.5))

  output_file.printf("JGLCeg = %.6g\n", soma.JGLCeg_glia_2013(0.5))

  output_file.printf("JGLCen = %.6g\n", soma.JGLCen_glia_2013(0.5))

  output_file.printf("JHKPFKn = %.6g\n", soma.JHKPFKn_glia_2013(0.5))

  output_file.printf("JHKPFKg = %.6g\n", soma.JHKPFKg_glia_2013(0.5))

  output_file.printf("ADPn = %.6g\n", soma.ADPn_glia_2013(0.5))

  output_file.printf("ADPg = %.6g\n", soma.ADPg_glia_2013(0.5))

  output_file.printf("JPGKn = %.6g\n", soma.JPGKn_glia_2013(0.5))

  output_file.printf("JPGKg = %.6g\n", soma.JPGKg_glia_2013(0.5))

  output_file.printf("JPKn = %.6g\n", soma.JPKn_glia_2013(0.5))

  output_file.printf("JPKg = %.6g\n", soma.JPKg_glia_2013(0.5))

  output_file.printf("JLDHn = %.6g\n", soma.JLDHn_glia_2013(0.5))

  output_file.printf("JLDHg = %.6g\n", soma.JLDHg_glia_2013(0.5))

  output_file.printf("JLACec = %.6g\n", soma.JLACec_glia_2013(0.5))

  output_file.printf("JLACcg = %.6g\n", soma.JLACcg_glia_2013(0.5))

  output_file.printf("JLACeg = %.6g\n", soma.JLACeg_glia_2013(0.5))

  output_file.printf("JLACen = %.6g\n", soma.JLACen_glia_2013(0.5))

  output_file.printf("Jmitoing = %.6g\n", soma.Jmitoing_glia_2013(0.5))

  output_file.printf("Jmitoinn = %.6g\n", soma.Jmitoinn_glia_2013(0.5))

  output_file.printf("Jmitooutn = %.6g\n", soma.Jmitooutn_glia_2013(0.5))

  output_file.printf("Jmitooutg = %.6g\n", soma.Jmitooutg_glia_2013(0.5))
  output_file.printf("Jshuttlen = %.6g\n", soma.Jshuttlen_glia_2013(0.5))

  output_file.printf("Jshuttleg = %.6g\n", soma.Jshuttleg_glia_2013(0.5))

  output_file.printf("JCKn = %.6g\n", soma.JCKn_glia_2013(0.5))

  output_file.printf("JCKg = %.6g\n", soma.JCKg_glia_2013(0.5))

  output_file.printf("JO2mcn = %.6g\n", soma.JO2mcn_glia_2013(0.5))

  output_file.printf("JO2mcg = %.6g\n", soma.JO2mcg_glia_2013(0.5))

  output_file.printf("JO2c = %.6g\n", soma.JO2c_glia_2013(0.5))

  output_file.printf("JGLCc = %.6g\n", soma.JGLCc_glia_2013(0.5))

  output_file.printf("JLACc = %.6g\n", soma.JLACc_glia_2013(0.5))

  output_file.printf("O2meanc = %.6g\n", soma.O2meanc_glia_2013(0.5))

  output_file.printf("hh_EL = %.6g\n", soma.hh_EL_glia_2013(0.5))

  output_file.printf("IL = %.6g\n", soma.IL_glia_2013(0.5))

  output_file.printf("hh_alpham = %.6g\n", soma.hh_alpham_glia_2013(0.5))

  output_file.printf("hh_betam = %.6g\n", soma.hh_betam_glia_2013(0.5))

  output_file.printf("hh_alphah = %.6g\n", soma.hh_alphah_glia_2013(0.5))

  output_file.printf("hh_betah = %.6g\n", soma.hh_betah_glia_2013(0.5))

  output_file.printf("hh_alphan = %.6g\n", soma.hh_alphan_glia_2013(0.5))

  output_file.printf("hh_betan = %.6g\n", soma.hh_betan_glia_2013(0.5))

  output_file.printf("hh_mCa = %.6g\n", soma.hh_mCa_glia_2013(0.5))

  output_file.printf("hh_tauh = %.6g\n", soma.hh_tauh_glia_2013(0.5))

  output_file.printf("hh_minfinity = %.6g\n", soma.hh_minfinity_glia_2013(0.5))

  output_file.printf("hh_ninfinity = %.6g\n", soma.hh_ninfinity_glia_2013(0.5))

  output_file.printf("hh_hinfinity = %.6g\n", soma.hh_hinfinity_glia_2013(0.5))

  output_file.printf("hh_taun = %.6g\n", soma.hh_taun_glia_2013(0.5))

  output_file.printf("INa = %.6g\n", soma.INa_glia_2013(0.5))

  output_file.printf("IK = %.6g\n", soma.IK_glia_2013(0.5))

  output_file.printf("ICa = %.6g\n", soma.ICa_glia_2013(0.5))

  output_file.printf("ImAHP = %.6g\n", soma.ImAHP_glia_2013(0.5))

  output_file.printf("Ipump = %.6g\n", soma.Ipump_glia_2013(0.5))

  output_file.printf("Fout = %.6g\n", soma.Fout_glia_2013(0.5))

  output_file.printf("ren = %.6g\n", soma.ren_glia_2013(0.5))

  output_file.printf("reg = %.6g\n", soma.reg_glia_2013(0.5))

  output_file.printf("rce = %.6g\n", soma.rce_glia_2013(0.5))

  output_file.printf("rcn = %.6g\n", soma.rcn_glia_2013(0.5))

  output_file.printf("rcg = %.6g\n", soma.rcg_glia_2013(0.5))

  output_file.printf("JGLCec = %.6g\n", soma.JGLCec_glia_2013(0.5))

  output_file.printf("JGLCgc = %.6g\n", soma.JGLCgc_glia_2013(0.5))

  output_file.printf("JGLCge = %.6g\n", soma.JGLCge_glia_2013(0.5))

  output_file.printf("JGLCne = %.6g\n", soma.JGLCne_glia_2013(0.5))

  output_file.printf("JLACce= %.6g\n", soma.JLACce_glia_2013(0.5))

  output_file.printf("JLACgc= %.6g\n", soma.JLACgc_glia_2013(0.5))

  output_file.printf("JLACge= %.6g\n", soma.JLACge_glia_2013(0.5))
  output_file.printf("JLACne= %.6g\n", soma.JLACne_glia_2013(0.5))

  output_file.printf("un= %.6g\n", soma.un_glia_2013(0.5))
  output_file.printf("ug= %.6g\n", soma.ug_glia_2013(0.5))

  output_file.printf("dAMPdATPn= %.6g\n", soma.dAMPdATPn_glia_2013(0.5))

  output_file.printf("dAMPdATPg= %.6g\n", soma.dAMPdATPg_glia_2013(0.5))

//  output_file.printf("vnstim= %.6g\n",soma.vnstim_glia_2013(0.5))

//  output_file.printf("vgstim= %.6g\n",soma.vgstim_glia_2013(0.5))


output_file.close()