Ortner_Simulation - SimTest.hoc

//Creation of SNc DA single compartment neuron with a PRE HH neuron activating an NMDA synapse 
// written by Antonios Dougalis, Dec 2014, Kozani
// extended by Antonios Dougalis, Aug 2015, Innsbruck
// finalised by Antonios Dougalis, Feb 2016, Ulm
//----------------------------------------------------------------------------
//  create DA compartment 
//----------------------------------------------------------------------------
load_file("nrngui.hoc")

create PRE 
PRE {
  diam=10
  L=10
  Ra = 35.4
  cm=1
  insert hh
}
  

create DA 
DA {
  diam=30
  L=30
  Ra = 35.4
  cm=1
  insert HHHmod
  insert A
  insert H
  insert gkca
  insert cabalanmod
  insert capumpmod
  insert calcha
  insert calchaT
  insert cal2
  insert cal3s
  insert cal3l
}

//----------------------------------------------------------------------------
//  insert IClamp in DA  terminal 
//----------------------------------------------------------------------------
access DA
objref stim      
stim = new IClamp(0.5)
stim.del = 0
stim.dur = 10000        
stim.amp = 0

objref stim1      
stim1 = new IClamp(0.5)
stim1.del = 2000
stim1.dur = 2000        
stim1.amp = 0

//----------------------------------------------------------------------------
//  insert IClamp in PRE terminal 
//----------------------------------------------------------------------------
access PRE
objectvar stim2      
PRE stim2 = new IClamp(0.5)
stim2.del = 4000
stim2.dur = 100        
stim2.amp = 0.2

//----------------------------------------------------------------------------
//  insert DA NMDA postsynaptic mechanism
//----------------------------------------------------------------------------
access DA
objectvar synNMDA
synNMDA = new NMDA(0.5)          // create synapse
DA synNMDA.loc(0.5)         // assign postsynaptic compartment
setpointer synNMDA.PRENMDA, PRE.v(0.5)    // assign presynaptic compartment

//----------------------------------------------------------------------------
//  insert DA AMPA and GABA postsynaptic mechanisms
//----------------------------------------------------------------------------

objectvar syn, syns
syn= new AMPA(0.5)
syns=new GABA(0.7)

//for AMPA receptor
objref synstim, stimA, trignetcon, netsynstim   //syn is a synapse, synstim is a NetStim object, stim is an IClamp object,
//triggernetcon & netsynstim are NetCon objects

synstim = new NetStims(0.5)       //use the name of the point process mod file that is NetStims
synstim.start = 0               //define its properties
synstim.interval = 1000
synstim.number = 100
synstim.noise = 1

stimA = new IClamp(0.5)           //define an IClamp that will be used to trigger the NetStims
stimA.del = 0
stimA.dur = 10000
stimA.amp = 0.0004

trignetcon = new NetCon(&IClamp[0.5].i,synstim,0.0003,0,1) //triger Netstim via IClamp using a NetCon
netsynstim = new NetCon(synstim,syn,0,0,0.00025) //make netsynstim a NetCon that connects the source NetStim[synstim] to the target syn

//for GABA receptor
objref synstims, stimG, trignetcons, netsynstims   //syn is a synapse, synstim is a NetStim object, stim is an IClamp object,
//triggernetcon & netsynstim are NetCon objects

synstims = new NetStimss(0.7)       //use the name of the point process mod file that is NetStimss
synstims.start = 0               //define its properties
synstims.interval = 500
synstims.number = 100
synstims.noise = 1

stimG = new IClamp(0.7)           //define an IClamp that will be used to trigger the NetStimss
stimG.del = 0
stimG.dur = 10000
stimG.amp = 0.0004

trignetcons = new NetCon(&IClamp[0.7].i,synstims,0.0003,0,1) //triger Netstim via IClamp using a NetCon
netsynstims = new NetCon(synstims,syns,0,0,0.0005) //make netsynstim a NetCon that connects the source NetStim[synstim] to the target syn



//----------------------------------------------------------------------------
//vbox utilities for graphing variables
//----------------------------------------------------------------------------
objref vbox, g, h, m
vbox = new VBox()
vbox.intercept(1)   //all following creations go into the "vbox" box
g = new Graph()
g.label("mV")
g.size(0,10000,-80,20)                                                                //g.size(xstart, xstop, ystart, ystop)
g.addvar("SNc DA NEURON voltage (mV)","DA.v(0.5)", 1,0)                            //g.addvar("label", "variable", color_index, brush_index)
g.flush()
graphList[0].append(g)
h = new Graph()
h.label("mV")
h.size(4000,4100,-80,20)                                                                //g.size(xstart, xstop, ystart, ystop)
h.addvar("PRESYNAPTIC HH NEURON voltage (mV)","PRE.v(0.5)", 2,0)                            //g.addvar("label", "variable", color_index, brush_index)
h.flush()
graphList[0].append(h)
m = new Graph()
m.label("mV")
m.size(4000,4200,-80,20)                                                                //g.size(xstart, xstop, ystart, ystop)
m.addvar("SNc DA NEURON voltage (mV)","DA.v(0.5)", 1,0)                            //g.addvar("label", "variable", color_index, brush_index)
m.flush()
graphList[0].append(m)
xpanel("")
xpanel()
vbox.intercept(0)   //ends intercept mode
vbox.map("Voltage Time Series",300, 100, -1, 1)          //draw the box and its contents and labels it and specify coordinates map("label", left, top, width, height)

objref vbox1, k, l
vbox1 = new VBox()
vbox1.intercept(1)   //all following creations go into the "vbox" box
k = new Graph()
k.label("mM")
k.size(0,10000,0,0.001)                                                                //g.size(xstart, xstop, ystart, ystop)
k.addvar("Intracellular Calcium (mM)","DA.cai(0.5)", 3,0)                            //g.addvar("label", "variable", color_index, brush_index)
k.flush()
graphList[0].append(k)
l = new Graph()
l.label("nA")
l.size(0,10000,-3,0)                                                                //g.size(xstart, xstop, ystart, ystop)
l.addvar("NMDA Current(nA)","NMDA[0].ica", 4,0)                            //g.addvar("label", "variable", color_index, brush_index)
l.flush()
graphList[0].append(l)
xpanel("")
xpanel()
vbox1.intercept(0)   //ends intercept mode
vbox1.map("Calcium and NMDA receptor",750, 100, -1, 1)          //draw the box and its contents and labels it and specify coordinates map("label", left, top, width, height)



objref vbox2, p, q
vbox2 = new VBox()
vbox2.intercept(1)   //all following creations go into the "vbox" box
p = new Graph()
p.label("mA/cm2")
p.size(3700,4200,0,0.1)                                                                //g.size(xstart, xstop, ystart, ystop)

                       //g.addvar("variable", color_index, brush_index)
p.addvar("Kdr","DA.ik_HHHmod(0.5)", 2,0)
p.addvar("A","DA.ik_A(0.5)", 3,0)
p.addvar("SK", "DA.ik_gkca(0.5)", 1, 0)
p.flush()
graphList[0].append(p)
q = new Graph()
q.label("mA/cm2")
q.size(3000, 4000,0,0.1)
q.addvar("DA all Calcium currents","DA.ica( 0.5 )", 1, 0)
q.addvar("T current","DA.ica_calchaT( 0.5 )", 2, 0)
q.addvar("N current","DA.ican_calcha( 0.5 )", 3, 0)
q.addvar("R current","DA.icar_calcha( 0.5 )", 4, 0)
q.addvar("Pump current","DA.icapump_capumpmod( 0.5 )", 5, 0)
q.addvar("All calcium again current","DA.ica_calcha( 0.5 )", 6, 0)
q.flush()
graphList[0].append(q)
xpanel("")
xpanel()
vbox2.intercept(0)   //ends intercept mode
vbox2.map("Pacemaking conductances",750, 100, -1, 1)          //draw the box and its contents and labels it and specify coordinates map("label", left, top, width, height)


//------------------------------------------------------------------------------------------------
// Management of Graphical user interface
//------------------------------------------------------------------------------------------------

{load_file("nrngui.hoc")}
objectvar save_window_, rvp_
objectvar scene_vector_[7]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}
{pwman_place(0,0,0)}

//Begin PointProcessGroupManager
{
load_file("pointgrp.hoc")
}
{
ocbox_ = new PointProcessGroupManager(0)
}
{object_push(ocbox_)}
{
}
{object_pop() doNotify()}
{
ocbox_ = ocbox_.v1
ocbox_.map("Electrode, Synapse & Noise Controller", 285, 750, 973.8, 396)
}
objref ocbox_
//End PointProcessGroupManager

{
xpanel("RunControl", 0)
v_init = -65
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 10000
xvalue("t","t", 2 )
tstop = 10000
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.025
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 8.95
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(0,115)
}
{
xpanel("SNc DA Neuron Parameter Box", 0)
xlabel("DA basal characteristics")
DA.L = 30
xvalue("L","DA.L", 1,"", 0, 0 )
DA.diam = 30
xvalue("diam","DA.diam", 1,"", 0, 0 )
DA.Ra = 35.4
xvalue("Ra","DA.Ra", 1,"DA.Ra += 0", 0, 1 )
DA.cm = 1
xvalue("cm","DA.cm", 1,"", 0, 0 )
DA.nseg = 1
xvalue("nseg","DA.nseg", 1,"", 0, 0 )
xlabel("")
xlabel("IONIC CONDUCTANCES")
DA.gnabar_HHHmod = 0.05
xvalue("gnabar_HHHmod","DA.gnabar_HHHmod", 1,"", 0, 0 )
DA.gkbar_HHHmod = 0.01
xvalue("gkbar_HHHmod","DA.gkbar_HHHmod", 1,"", 0, 0 )
DA.glbar_HHHmod = 8e-005
xvalue("glbar_HHHmod","DA.glbar_HHHmod", 1,"", 0, 0 )
DA.gkAbar_A = 0.01
xvalue("gkAbar_A","DA.gkAbar_A", 1,"", 0, 0 )
DA.ghbar_H = 0.0002
xvalue("ghbar_H","DA.ghbar_H", 1,"", 0, 0 )
DA.gkcabar_gkca = 0.0004
xvalue("gkcabar_gkca","DA.gkcabar_gkca", 1,"", 0, 0 )
DA.gcanbar_calcha = 5e-005
xvalue("gcanbar_calcha","DA.gcanbar_calcha", 1,"", 0, 0 )
DA.gcalbar_calcha = 0
xvalue("gcalbar_calcha","DA.gcalbar_calcha", 1,"", 0, 0 )
DA.gcarbar_calcha = 5e-005
xvalue("gcarbar_calcha","DA.gcarbar_calcha", 1,"", 0, 0 )
DA.gcatbar_calchaT = 7.5e-005
xvalue("gcatbar_calchaT","DA.gcatbar_calchaT", 1,"", 0, 0 )
DA.gcal2bar_cal2 = 0
xvalue("gcal2bar_cal2","DA.gcal2bar_cal2", 1,"", 0, 0 )
DA.gcal3sbar_cal3s = 0
xvalue("gcal3sbar_cal3s","DA.gcal3sbar_cal3s", 1,"", 0, 0 )
DA.gcal3lbar_cal3l = 0
xvalue("gcal3lbar_cal3l","DA.gcal3lbar_cal3l", 1,"", 0, 0 )
xlabel("")
xlabel("Ionic reversal potentials")
DA.el_HHHmod = -32
xvalue("el_HHHmod","DA.el_HHHmod", 1,"", 0, 0 )
DA.ena = 50
xvalue("ena","DA.ena", 1,"", 0, 0 )
DA.eh_H = -40
xvalue("eh_H","DA.eh_H", 1,"", 0, 0 )
DA.ek = -73
xvalue("ek","DA.ek", 1,"", 0, 0 )
xlabel("")
xpanel(1572,96)
}

xpanel("SNc DA Neuron Calcium Parameterisation Box", 0)
xlabel("")
xlabel("Calcium Parameterisation")
DA.kmkca_gkca = 0.00019
xvalue("kmkca_gkca","DA.kmkca_gkca", 1,"", 0, 0 )
DA.n_gkca = -4e-005
xvalue("n_gkca","DA.n_gkca", 1,"", 0, 0 )
DA.cainit_cabalanmod = 2e-005
xvalue("cainit_cabalanmod","DA.cainit_cabalanmod", 1,"", 0, 0 )
DA.fCa_cabalanmod = 0.025
xvalue("fCa_cabalanmod","DA.fCa_cabalanmod", 1,"", 0, 0 )
DA.icapumpmax_capumpmod = 0.01056
xvalue("icapumpmax_capumpmod","DA.icapumpmax_capumpmod", 1,"", 0, 0 )
DA.kmp_capumpmod = 0.0001
xvalue("kmp_capumpmod","DA.kmp_capumpmod", 1,"", 0, 0 )
DA.kmn_calcha = 0.0001
xvalue("kmn_calcha","DA.kmn_calcha", 1,"", 0, 0 )
DA.kml_calcha = 0.00045
xvalue("kml_calcha","DA.kml_calcha", 1,"", 0, 0 )
xpanel(1172,96)
objectvar scene_vector_[1]
{doNotify()}



//---------------------------------------------------------------------------------------------
// initiation of messaging
//---------------------------------------------------------------------------------------------
//boolean_dialog("Welcome to my SNc DA neuron model. Would you like to start simulating?", "YES", "NO")
//continue_dialog("Please read the Help File first before running simulations")

objref vbox3, vbox4
vbox3 = new VBox()
vbox3.intercept(1)   //all following creations go into the "vbox" box
xpanel("Welcome to my SNc DA neuron model", 1)
xlabel("")
xlabel("Would you like to start simulating?")
xlabel("")
xbutton("YES", "YES()")
xbutton("NO", "NO()")
xpanel()
vbox3.intercept(0)   //ends intercept mode
vbox3.map("Welcome to my SNc DA neuron model",600, 600, -1, 1)

proc vbox3UNMAP(){
vbox3.unmap
}

proc vbox4map(){
vbox4 = new VBox()
vbox4.intercept(1)   //all following creations go into the "vbox" box
//xpanel("Please read the Help File first before running simulations", 1)
//xlabel("Please read the Help File first before running simulations")
//xpanel(100,100)
continue_dialog("Please read the Help File first before running simulations")
//vbox4.dialog("Disclaimer", "Acc")
vbox4.intercept(0)   //ends intercept mode
//vbox4.dismiss_action("vbox4.unmap()")
//vbox4.map("test")
}

proc YES(){
init()
vbox3UNMAP()
vbox4map()
mkHELP()

}

proc NO(){
quit()
}


//----------------------------------------------------------------------------------------------
//creation and loading of a help file
//----------------------------------------------------------------------------------------------
proc mkHELP() {
xpanel("HELP_FILE")
xlabel("")
xlabel("Welcome to my SNc DA neuron model!") 
xlabel("-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------")
xlabel("You can minimise this help window at any time and find it again through the NEURON Main menu by clicking on the Window command")
xlabel("If you accidentally close the Help window, Click on File command of the NEURON Main Menu and choose load session then select")
xlabel("Help File.ses from your default folder where all simulation files are stored and it will reappear.")
xlabel("-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------")
xlabel("Brief Model description")
xlabel("The model consists of two neurons, the SNc DA neuron, named DA, which is endowed with a number of conductances and a Hogdkin and Huxley neuron,")
xlabel("named PRE, that controls the activation of an NMDA receptor which is located postsynaptically on the DA neuron.")
xlabel("The PRE neuron does not fire spontaneously and in order to activate the NMDA receptor located on DA we need to excite PRE via brief current clamp injections controlled by electrodes")
xlabel("-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------")
xlabel("How to introduce current clamp electrodes in the model")
xlabel("On the PointProcessGroupManager panel window, press and hold the Show All button and then select option IClamp.")
xlabel("Three electrodes with controllers for their parameters will automatically appear in that window.")
xlabel("They are labelled clearly to which entity they belong. Two are located on the DA neuron and one on the PRE neuron.")
xlabel("Use those to inject positive/negative current into the DA or the PRE neuron by selecting the ON time of current injection as a delay from the start of simulation, its duration and its amplitude.")
xlabel("-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------")
xlabel("How to introduce the NMDA receptor in the model")
xlabel("On the PointProcessGroupManager window, press and hold the Show All button and select the option NMDA.")
xlabel("The NMDA receptor, which is located on the DA neuron,will automatically appear on the window")
xlabel("NB: the IClamp on the PRE neuron is controlling the timing of the activation of the NMDA receptor. Use this electrode to activate the NMDA receptor")
xlabel("-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------")
xlabel("How to START running simulations")
xlabel("Use the RunControl box and press Init&Run. You should be able to observe the results of the simulations in the appropriate Graphic windows")
xlabel("The graphical windows are labelled with the parameter that is being monitored, i.e. DA.v(0.5) is DA neuron voltage.")
xlabel("DA.cai is the intracellular calcium and NMDA.ica is the current flowing through the NMDA receptor.")
xlabel("-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------")
xlabel("How to SAVE simulations")
xlabel("To export any of the resultant traces simply right click and hold on the graphical window's top left corner box and select Pick Vector")
xlabel("Then click and hold on the actual trace on the graphical window and the trace would go red")
xlabel("In the NEURON main Menu, select command Vector and click on Save to File. You will be prompt to give the file a name")
xlabel("Please name the file and give it an extension, name.dat. You have just exported the simulation trace as a dat file which will be saved in the default simulation folder")
xlabel("-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------")
xlabel("This model and files will be updated when necessary by the modeller and implementer.")
xlabel("For any queries or problems with this model please contact me on antonios.dougalis@uni-ulm.de.")
xlabel("When contacting me about a problem please make sure you include a short description/message AND a screenshot of any messages that you get in the command prompt and state briefly your actions.")
xlabel("This would help me enormously in addressing the problem in a timely manner.")
xlabel("")
xlabel("")
xlabel("Have fun !")
xlabel("")
xlabel("Dr. Antonios Dougalis")
xlabel("Molecular Neurophysiology")
xlabel("University of Ulm")
xlabel("")
xlabel("February 2016, Ulm")
xlabel("")
xpanel(200,100)
}

//boolean_dialog("Welcome to my SNc DA neuron model. Would you like to start simulating?", "YES", "NO") // or use it in this format coupled to the vbox b =  vbox.dialog("label", "Accept", "Cancel") 


//--------------------------------------------------------------------------------------------------------------------------------------------
//creating panels and functionalities for saving files
//--------------------------------------------------------------------------------------------------------------------------------------------

objref vbox5
vbox5 = new VBox()
vbox5.intercept(1)   //all following creations go into the "vbox" box
xpanel("Saving Simulation Data", 1)
xlabel("")
xlabel("Would you like to save simulated data?")
xlabel("")
xbutton("YES, Save my data", "savalldata()")
xbutton("Do Not Save Yet", "run()")
xbutton("Close dialogue, Do Not want to Save anything", "vbox5UNMAP()")

xpanel()
vbox5.intercept(0)   //ends intercept mode
vbox5.map("Saving Simulations",300, 0, -1, 1)

proc vbox5UNMAP(){
vbox5.unmap
}

objref rect, recv, recik, recica        //create vectors to hold t & v values from simulations
rect = new Vector(40001)
recv = new Vector(40001)
recik = new Vector(40001)
recica = new Vector(40001)          
rect.record(&t)                       //record t values in the vector rect 
recv.record(&DA.v(0.5))               //record v values in the vector recv
recik.record(&DA.ik(0.5))             //record ik values in the vector recik 
recica.record(&DA.ica(0.5))             //record ica values in the vector recica 

objref savdata              //create a file to be opened and saved and to hold values from vectors rect and recv
savdata = new File()
savdata.wopen("DAanalysis.dat")
//savdata.seek()
savdata.printf("t DA.v(0.5)\n")             //first line of data file
savdata.printf("ms mV mA/cm2 mA/cm2\n")       //second line of data for units
savdata.printf("time voltage ik ica\n")       //third line of data for labels      
//savdata.printf("%d\n",rect.size())
//for i=0, rect.size()-1 {
//    savdata.printf("%g\ %g\n", rect.x(i), recv.x(i), recik.x(i), recica.x(i))
//}


objref tempmatrix
tempmatrix = new Matrix(40001, 4)
tempmatrix.setcol(0, rect)
tempmatrix.setcol(1, recv)
tempmatrix.setcol(2, recik)
tempmatrix.setcol(3, recica)
tempmatrix.fprint(savdata, " %g")

proc savalldata(){
savdata.wopen("DAanalysis.dat")
//savdata.seek()
savdata.printf("t DA.v(0.5)\n")             //first line of data file
savdata.printf("s mV mA/cm2 mA/cm2\n")       //second line of data for units
savdata.printf("time voltage ik ica\n")       //third line of data for labels  
savdata.printf("%d\n",rect.size())
for i=0, rect.size()-1 {
   savdata.printf("%g\ %g\ %g\ %g\n", rect.x(i), recv.x(i), recik.x(i), recica.x(i))
}
tempmatrix.fprint(savdata, " %g")
savdata.close()
}

//--------------------------------------------------------------------------------------------------------------------
//creation of an extra Hbox containing a graph and a popup asking for further parameters to graphed be saved in files
//----------------------------------------------------------------------------------------------------------------------


objref graphix
objref choose

proc newgrapher(){
graphix = new Graph(0)
choose = new SymChooser()
choose.run()
}

objref Hbox1
Hbox1 = new HBox()
Hbox1.intercept(1)
xpanel ("")
xlabel("specify more variables to file")
xbutton ("YES","newgrapher()")
xbutton("NO", "NO()")
xpanel()
Hbox1.intercept(0)
Hbox1.map

proc NO(){
Hbox1.unmap()
}