// Created by FF (2018)
// Functions here are used to calculate extracellular conductances for PNs and ITNs
func calcconduc_PN() { local x_pos,y_pos,z_pos,sec,dmin,dmin2,H,r_sq,ds_comp,Ls,h,l,phi,conduct,r localobj x_end,x_comp,d_end,orilocal,locationlocal,eleclocal
x_end = new Vector(3)
d_end = new Vector(3)
x_comp = new Vector(3)
dmin = diam_soma_p/2 + diam_shank/2 + extralimit
dmin2=dmin*dmin ///set limit for electrodes
sec=$2 //sec=0 for soma,1-8 for adend, 9-15 for pdend
if (sec==0) {
x_pos=$o1.x[0] y_pos=$o1.x[1] z_pos=$o1.x[2] ///location of soma
xe=$o4.x[0] ye=$o4.x[1] ze=$o4.x[2] ///location of elec
dis_soma=$5
if (dis_soma<dmin) {
dis_soma=dmin
}
r = dis_soma
} else if (sec>=9) { ///calculate for pdend
ds_comp = abs(pdend_L_p/nseg_pdend_p)
x_end.x[0] = $o1.x[0]-((sec-nseg_adend_p)*ds_comp+diam_soma_p/2)*$o3.x[0]
x_end.x[1] = $o1.x[1]-((sec-nseg_adend_p)*ds_comp+diam_soma_p/2)*$o3.x[1]
x_end.x[2] = $o1.x[2]-((sec-nseg_adend_p)*ds_comp+diam_soma_p/2)*$o3.x[2]
x_comp.x[0] = $o3.x[0]*(-ds_comp)
x_comp.x[1] = $o3.x[1]*(-ds_comp)
x_comp.x[2] = $o3.x[2]*(-ds_comp)
d_end.x[0] = $o4.x[0]-x_end.x[0]
d_end.x[1] = $o4.x[1]-x_end.x[1]
d_end.x[2] = $o4.x[2]-x_end.x[2]
H = (d_end.x[0]*x_comp.x[0]+d_end.x[1]*x_comp.x[1]+d_end.x[2]*x_comp.x[2])/ds_comp
r_sq = (d_end.x[0]*d_end.x[0]+d_end.x[1]*d_end.x[1]+d_end.x[2]*d_end.x[2])-H*H
Ls = H + ds_comp
if (H*Ls>0) {
if (abs(H)>abs(Ls)) {
h = abs(H)
l = abs(Ls)
} else {
h = abs(Ls)
l = abs(H)
}
if (l*l+r_sq<dmin2) {
l = sqrt(dmin2-r_sq) //set limit
h = l+ds_comp
//n = n+1;
}
phi = log((sqrt((h*h)+r_sq)+h)/(sqrt((l*l)+r_sq)+l))
} else {
if (r_sq<dmin2) { //set limit
r_sq = dmin2
}
phi = log(((sqrt((Ls*Ls)+r_sq)+Ls)*(sqrt((H*H)+r_sq)-H))/r_sq)
}
r = ds_comp/phi
} else { ////calculate for a_dend
ds_comp = abs(adend_L_p/nseg_adend_p)
x_end.x[0] = $o1.x[0]+(sec*ds_comp+diam_soma_p/2)*$o3.x[0]
x_end.x[1] = $o1.x[1]+(sec*ds_comp+diam_soma_p/2)*$o3.x[1]
x_end.x[2] = $o1.x[2]+(sec*ds_comp+diam_soma_p/2)*$o3.x[2]
x_comp.x[0] = $o3.x[0]*(ds_comp)
x_comp.x[1] = $o3.x[1]*(ds_comp)
x_comp.x[2] = $o3.x[2]*(ds_comp)
d_end.x[0] = $o4.x[0]-x_end.x[0]
d_end.x[1] = $o4.x[1]-x_end.x[1]
d_end.x[2] = $o4.x[2]-x_end.x[2]
H = (d_end.x[0]*x_comp.x[0]+d_end.x[1]*x_comp.x[1]+d_end.x[2]*x_comp.x[2])/ds_comp
r_sq = (d_end.x[0]*d_end.x[0]+d_end.x[1]*d_end.x[1]+d_end.x[2]*d_end.x[2])-H*H
Ls = H + ds_comp
if (H*Ls>0) {
if (abs(H)>abs(Ls)) {
h = abs(H)
l = abs(Ls)
} else {
h = abs(Ls)
l = abs(H)
}
if (l*l+r_sq<dmin2) {
l = sqrt(dmin2-r_sq) //set limit
h = l+ds_comp
//n = n+1;
}
phi = log((sqrt((h*h)+r_sq)+h)/(sqrt((l*l)+r_sq)+l))
} else {
if (r_sq<dmin2) { //set limit
r_sq = dmin2
}
phi = log( ((sqrt((Ls*Ls)+r_sq)+Ls)*(sqrt((H*H)+r_sq)-H))/r_sq )
}
r = ds_comp/phi
}
conduct = 1/(4*PI*r*sigma)*0.1 //LFP unit is 100uV.
return conduct
}
func calcconduc_ITN() { local x_pos,y_pos,z_pos,sec,dmin,dmin2,H,r_sq,ds_comp,Ls,h,l,phi,conduct localobj x_end,x_comp,d_end,orilocal,locationlocal,eleclocal
x_end = new Vector(3)
d_end = new Vector(3)
x_comp = new Vector(3)
dmin = diam_soma_I/2 + diam_shank/2 + extralimit dmin2=dmin*dmin ///set limit for electrodes
sec=$2 //sec=0 for soma,otherwise for adend
if (sec==0) {
x_pos=$o1.x[0] y_pos=$o1.x[1] z_pos=$o1.x[2] ///location of soma
xe=$o4.x[0] ye=$o4.x[1] ze=$o4.x[2] ///location of elec
dis_soma=$5
if (dis_soma<dmin) {
dis_soma=dmin
}
r = dis_soma
} else { ////calculate for dend
ds_comp = abs(dend_L_I/nseg_dend_I)
x_end.x[0] = $o1.x[0]+(sec*ds_comp+diam_soma_I/2)*$o3.x[0]
x_end.x[1] = $o1.x[1]+(sec*ds_comp+diam_soma_I/2)*$o3.x[1]
x_end.x[2] = $o1.x[2]+(sec*ds_comp+diam_soma_I/2)*$o3.x[2]
x_comp.x[0] = $o3.x[0]*(ds_comp)
x_comp.x[1] = $o3.x[1]*(ds_comp)
x_comp.x[2] = $o3.x[2]*(ds_comp)
d_end.x[0] = $o4.x[0]-x_end.x[0]
d_end.x[1] = $o4.x[1]-x_end.x[1]
d_end.x[2] = $o4.x[2]-x_end.x[2]
H = (d_end.x[0]*x_comp.x[0]+d_end.x[1]*x_comp.x[1]+d_end.x[2]*x_comp.x[2])/ds_comp
r_sq = (d_end.x[0]*d_end.x[0]+d_end.x[1]*d_end.x[1]+d_end.x[2]*d_end.x[2])-H*H
Ls = H + ds_comp
if (H*Ls>0) {
if (abs(H)>abs(Ls)) {
h = abs(H)
l = abs(Ls)
} else {
h = abs(Ls)
l = abs(H)
}
if (l*l+r_sq<dmin2) {
l = sqrt(dmin2-r_sq) //set limit
h = l+ds_comp
}
phi = log((sqrt((h*h)+r_sq)+h)/(sqrt((l*l)+r_sq)+l))
} else {
if (r_sq<dmin2) { //set limit
r_sq = dmin2
}
phi = log( ((sqrt((Ls*Ls)+r_sq)+Ls)*(sqrt((H*H)+r_sq)-H))/r_sq )
}
r = ds_comp/phi
}
conduct = 1/(4*PI*r*sigma)*0.1 //LFP unit is 100uV.
return conduct
}
///////////////////////////////////////////////////////////