/* HVCRA_neuron.hoc
Model for HVC(RA) projection neuron with an unmyelinated axon.
Simple ball and stick model of axon and soma. The axon
has two types of sections one contains a mitodondrion and the other does not.
The section with the mitochondrion has increased intracellular resisitivity.
Add sodium and kv1 currents to axon from model of Cohen et al 2020
Keep soma passive. Increase channel densities in axon.
Modified to run in batch mode and compute conduction velocity.
Run for 40 deg C and axon diameters of 0.2, 0.4, and 0.6 um.
Model with additional 9 micron initial mitochondrion free section.
Conduction velocity computed in the middle of the axon to avoid
sealed end effects.
The parameters that set the topology of the axon for a given longitudinal mitochondrial
coverage are specified in the hoc file that loads this file. This includes the
number of axon sections, the length and nseg of the mitochondria free sections and the
axon sections used to compute conduction velocity. The length of mitochondrion containing sections
are fixed at 1 micron. The axon diameter is selected in the gui, as is the longitudinal mitochondrial coverage.
This simulation computes conduction velocity for the cross-sectional mitochondrial occupancy
p ranging from 0.1 to 0.8. The calculation of conduction velocity for the mitochondria free axon, p=0,
is included for completeness, but was computed in a different simulation. The conduction
velocity is used calculate the additional time needed for the AP to travel 3 mm due to mitochondria in the
axon. The conduction velocity (m/s) and additional delay (ms) are plotted for the cross-sectional mitochondrial
occupancy, p.
*/
//parameters
proc model_globels() {
celsius = 40 //for birds temperature 38-43 deg C
p=$1 // proportion of crosssection occupied by mitochondrion
rhoax = 100 // cytoplasmic resistivity Ohm cm
rhomit = 10000 // mitochondrion resistivity
somaL=6 // soma length = soma diameter
somaD=6
//axonD=0.4 axon diameter selected in the gui
//axonL= 7 length of axon section without mitochondrion depnds on longitudinal coverage
mitochdL = 1 // length of section with mitochondrion
gleakpas = .0003
eleakpas = -70
//cm = 1
}
/*
//topological parameters
nmito=250 // number of mitochondria sections
naxon=250 // number of axon sections
*/
// Build neuron topology
create soma, axonint, axon[naxon], mitochd[nmito]
connect axonint(0), soma(1)
connect mitochd[0](0), axonint(1)
for i=0, nmito-2 {
connect axon[i](0), mitochd[i](1)
connect mitochd[i+1](0), axon[i](1)
}
connect axon[naxon-1](0), mitochd[nmito-1](1)
// Add biophysical mechanisms
proc morph() {
p=$1 //cross-sectional mitochondrial occupancy
axonD=$2 //axon diameter
soma {
nseg=3
Ra=100
cm=1
L=somaL
diam=somaD
insert pas
g_pas = gleakpas
e_pas = eleakpas
}
axonint{
nseg=11
Ra=100
cm=1
L= 9 //(um) initial axon section
diam=axonD
insert pas
g_pas = 10*gleakpas
e_pas = eleakpas
insert na
gbar_na = 1000 //pS/cm2
insert kv1
gbar_kv1 = 3000 //ps/cm2
}
for i=0,naxon-1 { //set parameters for sections without mitochondria
axon[i] {
nseg=nseg_axon //nseg adjusted to keep numerical compartments < 1 micron
Ra=100
cm=1
L=axonL
diam=axonD
insert pas
g_pas = 10*gleakpas
e_pas = eleakpas
insert na
gbar_na = 1000 //pS/cm2
insert kv1
gbar_kv1 = 3000 //ps/cm2
}
}
for i=0,nmito-1 { //set parameters for sections with mitochondria
mitochd[i] {
nseg=3
Ra= rhoax*rhomit/(p*rhoax+(1-p)*rhomit)
cm=1
L=mitochdL
diam=axonD
insert pas
g_pas = 10*gleakpas
e_pas = eleakpas
insert na
gbar_na = 1000 //pS/cm2
insert kv1
gbar_kv1 = 3000 //ps/cm2
}
}
}
access soma
// current pulse stimulus
objref stim, v1, v2, tvec, g1, g2, b1
soma stim = new IClamp(0.5)
stim.amp = 0.5
stim.del = 10
stim.dur = 0.5
v1 = new Vector(20000)
v2 = new Vector(20000)
tvec = new Vector(20000)
tvec.record(&t)
// Set up arrays for plots of conduction velocity and additional delay
double pp[9], condvel[9], delay[9]
proc init1(){
finitialize(v_init)
fcurrent()
}
proc integrate1() {
while (t<tstop ) {
fadvance()
}
}
proc run1() {
p=$1
axonD=$2
tstop=50
dt=.0025
v_init=-70
model_globels(p, axonD)
morph(p, axonD)
v1.record(&axon[naxon1].v(0.5))
v2.record(&axon[naxon2].v(0.5))
init1(v_init)
integrate1(tstop)
//soma for (x) fprint("%.10g \n", v(x))
//neurite for (x) fprint("%.10g %.10g %.10g \n", v(x), m_hhfxd(x), h_hhfxd(x))
}
//interpolated time at which voltage=threshold
func where() {local i, v1, v2, t1, t2, x
i = $o1.indwhere(">=", $2)
v1 = $o1.x[i-1]
v2 = $o1.x[i]
t1 = tvec.x[i-1]
t2 = tvec.x[i]
x = ($2 - v1)/(v2-v1)
return t1+x*(t2-t1)
}
proc velocity() { local i, thresh // m/s
p=$1
axonD=$2
thresh = -5
run1(p, axonD)
t1 = where(v1, thresh)
t2 = where(v2, thresh)
delt=t2-t1
delx=200/1000 //mm length
vel = delx/delt //compute velocity
//fprint ("%.10g %.10g %.10g %.10g %.10g %.10g %.10g \n", axonD, p, t1, t2, delt, delx, vel)
//printf ("%.10g %.10g %.10g \n", axonD, p, vel)
}
strdef file
proc batch2(){local i, ii, x, z, v
// set up external output file to save conduction velocity
/* v=1
sprint(file,"HVCRA_outputv%d.dat",v)
wopen(file)
fprint("diam p time1 time2 delt delx velocity \n") //print headings
*/
printf ("%.10s %.10s %.10s %.10s\n", "diam", "p", "velocity", "delay") //print to screen
celsius=40
// Conduction velocity for mitochondria-free axon was computed in a separate simulation
// and depends on axon diameter. The conduction velocity for p=0 and axonD is set below.
i=0
pp[i] = 0 //p=0
if (axonD == 0.2) {
condvel[i]= 0.31228 //m/s
}
if (axonD == 0.4) {
condvel[i]= 0.44188 //m/s
}
if (axonD == 0.6) {
condvel[i]= 0.54132 //m/s
}
delay0 = 0.003/condvel[i] // time needed for AP to travel 3mm in mitochondrion-free axon
delay[i]= (0.003/condvel[i]-delay0)*1000 // additonal delay(ms) due to mitochondria
printf ("%.10g %.10g %.10g %.10g \n", axonD, pp[i], condvel[i], delay[i])
i=i+1
for (x=1.0; x<=8.0; x=x+1.0) {
p=0.1*x // proportion of cross-section occupied by mitochrondrion
pp[i]=p
velocity(p, axonD)
condvel[i]=vel
delay[i]= (0.003/condvel[i]-delay0)*1000 // additonal delay(ms) due to mitochondria
printf ("%.10g %.10g %.10g %.10g \n", axonD, pp[i], condvel[i], delay[i])
i=i+1
}
// Plot conduction velocity (m/s) vs mitochondrial occupancy, p
// and plot additional delay (ms) vs mitochondrial occupancy, p
b1 = new VBox()
b1.intercept(1)
g1 = new Graph()
g2 = new Graph()
g1.size(0,0.8,0.20,0.55)
g1.beginline()
for (ii=0; ii<=8; ii=ii+1) {
g1.line(pp[ii], condvel[ii])
}
g1.flush()
g2.size(0,0.8,0.0,3.5)
g2.beginline()
for (ii=0; ii<=8; ii=ii+1) {
g2.line(pp[ii], delay[ii])
}
g2.flush()
b1.intercept(0)
b1.map("Conduction velocity(m/s), delay(ms)")
//wopen() close output file
}
batch2()
/*proc timer() {
startsw()
batch2()
print "the elapsed time is ", stopsw(), "secs. \n"
}*/