I was wondering if you could assist my understanding when simulating a neuron with multiple synaptic inputs along its dendrite.
Let me describe my setup and my code to give the bigger picture, followed by the questions I have.
I am using Zilany's Model 2014 to generate a number of spike trains (i.e. you specify the number of ANFs = ANFx) for a number of center frequencies (i.e. CFx). Each generated spike train consists of an array of time instances, but it must be noted that these arrays are not of equal length/have the same number of elements per array.
The output of Zilany's Model is saved as such:
Code: Select all
st1,st2,st3,....,stn
________________________________
0 C | ANF1: [t1, t2, t3, ... , tn] |
1 F | ANF2: [t1, t2, t3, ... , tn] |
2 : | ANF3: [t1, t2, t3, ... , tn] |
3 1 | ANF4: [t1, t2, t3, ... , tn] |
4 | ANFx: [t1, t2, t3, ... , tn] |
________________________________
: C | ANF1: [t1, t2, t3, ... , tn] |
: F | ANF2: [t1, t2, t3, ... , tn] |
: : | ANF3: [t1, t2, t3, ... , tn] |
: 2 | ANF4: [t1, t2, t3, ... , tn] |
: | ANFx: [t1, t2, t3, ... , tn] |
________________________________
: C | ANF1: [t1, t2, t3, ... , tn] |
: F | ANF2: [t1, t2, t3, ... , tn] |
: : | ANF3: [t1, t2, t3, ... , tn] |
: x | ANF4: [t1, t2, t3, ... , tn] |
n | ANFx: [t1, t2, t3, ... , tn] |
________________________________
The reason for the explanation of the input to my NEURON model is to understand the setup of the neuron, more specifically the dendrite. The "nseg" value for the dendrite of my model is dependent on the number of centre frequencies (CF) that I choose in Zilany. The reason being is, I want each segment of the dendrite to be responsible for a single CF (i.e. I want to keep the tonotopic arrangement of the termination of the spike trains along the dendrite).
The NEURON model I created is that of a simple neuron (i.e. one soma, one dendrite, one axon) for testing purposes before I implement my own biophysically defined neuron. The code below shows that setup.
Code: Select all
"""Single Simple Neuron:"""
#Create Cell Sections:
soma = h.Section(name = 'soma') #This is the soma.
dend1 = h.Section(name = 'dend1') #This is the dendrite.
axon = h.Section(name = 'axon') #This is the axon
#Define Cell Topology:
dend1.connect(soma(0), 0) #The dendrite end 0 connects to the soma end 0.
axon.connect(soma(1),) #The axon end (default 0) connects to the soma end 1.
#Define Geometry of Sections of the Cell:
"""Dendrite1"""
dend1.diam = 3 #Diameter [um]
dend1.L = 250 #Length [um]
dend1.nseg = len(CF_Range) #Number of Segments equals the number of CFs used in Zilany. NB must always be ODD!
"""Soma"""
soma.diam = 25 #Diameter [um]
soma.L = 25 #Length [um]
soma.nseg = 1 #Number of Segments
"""Axon"""
axon.diam = 3 #Diameter [um]
axon.L = 32 #Length [um]
axon.nseg = 1 #Number of Segments
h.define_shape()
#Define the Biophysics of each Section of the Cell:
for sec in h.allsec():
sec.Ra = 100 # Axial resistance in Ohm * cm
sec.cm = 0.9 # Membrane capacitance in micro Farads / cm^2
"""Insert active Hodgkin-Huxley current in the soma:"""
soma.insert('hh') #Inserts HH biophysics in soma.
for seg in soma: #For all defined segments of the soma insert the following:
seg.hh.gnabar = 0.12 # Sodium conductance in S/cm2
seg.hh.gkbar = 0.036 # Potassium conductance in S/cm2
seg.hh.gl = 0.0003 # Leak conductance in S/cm2
seg.hh.el = -54.3 # Reversal potential in mV
"""Insert active Hodgkin-Huxley current in the dendrite:"""
dend1.insert('hh') #Inserts HH biophysics in dendrite
for seg in dend1: #For all defined segments of the dendrite insert the following:
seg.hh.gnabar = 0.12 # Sodium conductance in S/cm2
seg.hh.gkbar = 0.036 # Potassium conductance in S/cm2
seg.hh.gl = 0.0003 # Leak conductance in S/cm2
seg.hh.el = -54.3 # Reversal potential in mV
# seg.pas.g = 0.001 # Passive conductance in S/cm2
# seg.pas.e = -65 # Leak reversal potential mV
"""Insert passive current in the axon:"""
axon.insert('pas') #Inserts pas biophysics in axon
for seg in axon: # For all defined segments of the axon insert the following:
seg.pas.g = 0.001 # Passive conductance in S/cm2
seg.pas.e = -65 # Leak reversal potential mV
h.topology() # Shows topology of created neuron
ps = h.PlotShape() # False tells h.PlotShape not to use NEURON's gui. Shows/draws the generated Octopus Cell.
ps.show(0)
My calculation of the placement array is as follows:
- nseg = number of CFs chosen (ALWAYS ODD)
- L = Length of each segment = 1/nseg
- HLen = Half the length of each segment = 0.5 * (1/nseg)
- Cfperseg = Number of segments per CF = nseg/number of CFs chosen
- CFs chosen = [CF1, CF2, ..., CFn]
- Placement = HLen + (x - 1)*(L), where x is the index in the list called CFs chosen
First step that I took was to restructure the output from the Zilany Model by taking into consideration that all spike trains related to a specific center frequency, which in turn related to a specific place/segment along the dendrite. Secondly, I also considered that I may use more than one dendrite, thus I may have to split the number of spike trains per center frequency evenly among the number of dendrites I specify/create. At the moment the model uses only one dendrite, but this is a future consideration. Therefore I ensured that the number of ANFs I specify in the Zilany Model is ALWAYS a multiple of the number of dendrites I create in NEURON.
Based on the above considerations I restructured the Zilany spiketrains such that I group all time elements per array index together, if the index does not exist in the original Zilany spike train output (i.e. as mentioned above each spike train does not contain the same number of elements), I append a value of 0. Therefore my spike train data I use for my NEURON model looks like this (i.e. a python List of Lists of Lists):
Code: Select all
Place on Dendrite: | 1st nsegment on each dendrite = CF1 || 2nd nsegment on each dendrite = CF2 || nth nsegment on each dendrite = CFn |
Dendrite number: | | D1 | ... | nth D | || | D1 | ... | nth D | || | D1 | ... | nth D | |
0 Simrun 1 = st1 [ [ [t1, t1, ...... ,nth t1 ] ... [t1, t1, ..... , nth t1 ] ][ [t1, t1, ...... ,nth t1 ] ... [t1, t1, ..... , nth t1 ] ][ [t1, t1, ...... ,nth t1 ] ... [t1, t1, ..... , nth t1 ] ] ]
1 Simrun 2 = st2 [ [ [t2, t2, ...... ,nth t2 ] ... [t2, t2, ..... , nth t2 ] ][ [t2, t2, ...... ,nth t2 ] ... [t2, t2, ..... , nth t2 ] ][ [t2, t2, ...... ,nth t2 ] ... [t2, t2, ..... , nth t2 ] ] ]
: : : : : : : : : :
n Simrun n = stn [ [ [tn, tn, ...... ,nth tn ] ... [tn, tn, ..... , nth tn ] ][ [tn, tn, ...... ,nth tn ] ... [tn, tn, ..... , nth tn ] ][ [tn, tn, ...... ,nth tn ] ... [tn, tn, ..... , nth tn ] ] ]
This follows the logic I used for the setup of the NEURON simulation as follows:
- For each segment along the dendrite to which I intend to input a synaptic input, I create an h.Exp2Syn() and append that to an h.List(). Thus each segment along the dendrite would have its own unique h.Exp2Syn() point process.
- I also create one h.VecStim() for each time instance per segment along each dendrite. I also append each h.VecStim() to a h.List()
- For each h.VecStim() created I create an h.Vector() to which the specified spike time per dendrite per segment of dendrite is inserted and that h.VecStim() can play.
- I then create a h.NetCon() object for all the h.VecStim() using the h.Exp2Syn() for that specific segment of the dendrite.
- I then set up the recording vectors for the tvec and voltage membrane of the soma.
- Finally, for each simrun I call h.run() and then I reinitialise repeating steps 1 to 4 before h.run() again.
Code: Select all
"""NEURON Simulation Run Setup:"""
#++++++++++++++++++++++++++++++++++++++
#Simulation Environmental Constants:
#++++++++++++++++++++++++++++++++++++++
tau1 = 0.07 #DECAY TIME IN ms.
tau2 = 0.34 #RISE TIME IN ms.
e = 0.0 #REVERSAL POTENTIAL IN mV.
weight = 2 #NetCon Weight [nS].
delay = 0 #NetCon Delay = 0.
threshold = 0 #NetCon Threshold = 0.
Temp = 34 #TEMPERATURE OF ENVIRONMENT.
T_Step = 0.005 #TIME STEP.
v0 = -63 #INITIAL VOLTAGE MEMBRANE IN mV.
Tstop = 10 #TIME OF STIMULATION IN s.
V_init = -63 #Initial voltage in mV of membrane.
Vm_soma = [] #List to append Vm Soma results to per simrun.
tvec = [] #List to append time vec results to per simrun.
Syn_I = [] #List to append Synaptic Current results to per simrun.
ALL_synI = [] #List to append ALL synaptic current results to.
for simrun in simrunList: #Per simrun/ number of elements in simrunList (A list of lists) execute the following:
h.celsius = Temp #Sets the NEURON Environment temperature
h.dt = T_Step #Sets the NEURON stimulation time step
h.v_init = V_init #Sets the NEURON initial voltage
h.finitialize(v0) #Initializes the NEURON initial voltage
exp2synList = h.List() #Creates a h.list to append all exp2syn created.
vecstimList = h.List() #Creates a h.list to append all vecstim created.
netconList = h.List() #Creates a h.list to append all netcons created.
synIList = h.List() #Creates a h.list to append all synI created.
for placealongdend in simrun: #For each segment along the dendrite per simrun execute the following:
dend = 1
for dendnum in placealongdend: #For each dendrite at the specified segment per simrun execute the following:
if dend == 1: #If the dend num = dendrite 1 execute the following.
dendtocall = dend1(Placement[place]) #Create variable dendtocall which specifies the neuron section and segment
elif dend == 2: #else If the dend num = dendrite 2 execute the following.
dendtocall = dend2(Placement[place]) #Create variable dendtocall which specifies the neuron section and segment
elif dend == 3: #else If the dend num = dendrite 3 execute the following.
dendtocall = dend3(Placement[place]) #Create variable dendtocall which specifies the neuron section and segment
elif dend == 4: #else If the dend num = dendrite 4 execute the following.
dendtocall = dend4(Placement[place]) #Create variable dendtocall which specifies the neuron section and segment
else: #else execute the following.
print("Number of dendrites != num of dends specified! CHECK CODE")
pass
"""setup Exp2Syn""" #The EXP2SYN setup
synapse = h.Exp2Syn(dendtocall) #Create h.Exp2Syn() point process for specific dendrite at specific dendrite segment.
exp2synList.append(synapse) #Append synapse to exp2synList.
exp2synList[-1].tau1 = tau1 #Set decay time.
exp2synList[-1].tau2 = tau2 #Set rise time.
exp2synList[-1].e = e #Set reversal potential.
"""Setup Vecstim""" #The VECSTIM setup
for time in dendnum: #For each time element in list execute the following:
vecstim = h.VecStim() #Create h.VecStim() for time element.
vecstimList.append(vecstim) #Append vecstim to vecstimList.
spiketime = h.Vector([time]) #Create a h.Vector() and insert the time element in the vector.
vecstimList[-1].play(spiketime) #Play the spiketime.
"""Setup Netcon:""" #The NETCON setup
net = h.NetCon(vecstimList[-1], exp2synList[-1]) #Create h.NetCon() using same exp2syn for all time elements in List per dendrite per segment of dendrite.
netconList.append(net) #Append net to netconList.
netconList[-1].weight[0] = weight #Set Netcon weight.
netconList[-1].delay = delay #Set Netcon delay.
netconList[-1].threshold = threshold #Set Netcon threshold.
dend += 1 #Increment dendrite number.
"""Run Simulation"""
somaR = h.Vector() #Create a vector for soma results.
somaR.record(soma(0.5)._ref_v) #Record the vm of the soma.
timeR = h.Vector() #Create a vector for time results.
timeR.record(h._ref_t) #Record the simulation time.
for es in exp2synList: #For each exp2syn created execute the followiing:
synI = h.Vector() #Create a vector for the synaptic current results
synIList.append(synI) #Append the synaptic current vectors to synIList
synIList[-1].record(es._ref_i) #Record the synaptic current.
h.tstop = Tstop #Set the simulation duration.
h.run() #Run NEURON simulation.
Vm_soma.append(np.array(somaR)) #Append Soma results to Vm_soma List
tvec.append(np.array(timeR)) #Append time results to tvec List
for sI in synIList: #For all synapse results in synIList execute the following:
Syn_I.append(np.array(sI)) #Append all synaptic current results to Syn_I List
ALL_synI.append(Syn_I) #Append Syn_I to ALL_synI List.
print("simulation Done") #print to cmd simulation is done.
- Seeing that I am only interested in the interspike intervals of the soma response, should I be inserting a h.NetCon() to record the tvec at the soma?
- Should I be re-initialising the parameters per Simrun, as seen in the above code? I ask, because in my mind by reinitializing I make all the parameters in the simulation = 0 or their equivalent initial values, which would basically mean a new cell per run. Am I correct? As in real life all these spike times would be arriving at the same cell, however consecutive action potentials will produce if threshold at the soma is exceeded and this is also dependent on how fast/slow the soma's ability is to polarize and depolarize based on its ionic channels.
- Is there not a way to cut out the restructuring of the Zilany Model Output, i.e. h.Vector([each spike train per CF]) or am I correct in my setup and use of h.VecStim and h.Vector, seeing that timing characteristics are important? Trying to improve computational efficiency I guess.
- Is appending 0 when restructuring the Zilany model as I mentioned above not influencing my data as an Exp2Syn() would occur at time t = 0, when in actuality no synaptic input or arrival of spike time on the specific segment of the dendrite should occur? Thus, is there not a NEURON shorthand telling NetCon/VecStim not to produce an event if the spike time = 0? The only way I have thought about this is simply implementing an "IF" "ELSE" statement to basically skip over the creation of a h. VecStim() if t = 0.
- Does NetCon.weight not need to vary along the dendrite length? I guess I am asking cause I understand that the weight variable measured in nS is the weight or for better words the strength of the Exp2Syn point process. Thus the further away from the soma the more weight should be applied to produce an action potential of same strength to an action potential produce by a synaptic input closer to the soma. Or is this totally not necessary.
- Furthermore, eventhough I have read up on h.NetCon(), I do not fully understand Netcon.threshold and how this influences the synaptic input. A little bit of enlightenment would be appreciated?
- Lastly, I think my biggest concern is, IS the synaptic instance along the dendrite being applied as I expect it should (seeing I'm newish in the use of NEURON and definitely new to coding using h.List())? Namely, have I forgotten/ not forgotten to clear a pointer or perhaps forgotten to consider something that will potentially affect the production of action potentials at my Soma? I do believe, that my output of my created model to be correct based on my results when I do plot my synaptic currents and voltage membrane potential of soma. Furthermore, I have checked that for each NetCon() object I am using the correct segment, source and target, as well as checking that each Exp2Syn() point process that I applied has the specified tau1, tau2, e and what the current it would be.
Regards