CellSwc_Ver2.hoc

//------------------------------------
// fixed 2015.01.08
//------------------------------------



{load_file("stdlib.hoc")}
/*
Use ExpSid to make synapse
*/

begintemplate CellSwc
public init, shape, translation, subsets, geom, biophys, geom_nseg, biophys_inhomo, makeSyn, makeNetCon
public Scalex, ScaleY, ScaleZ, Shrinkagex, Shrinkagey, Shrinkagez
public all, secInput, secOutput
public SynList
public Dend, SectionNum

public synlist,synapse
public SynNumberList, numSynapse
public mkSynapses

public synlist_,mkSyn

public mkIclamp
// will be re-created in shape()

public ID, Type, pID //pID -> parent ID

create Dend[1]

//------------------------------------------------------------------------------
// swc header informations
//------------------------------------------------------------------------------
strdef HeaderOriginalSource
strdef HeaderCreature
strdef HeaderRegion
strdef HeaderFieldLayer
strdef HeaderType
strdef HeaderContributor
strdef HeaderReference
strdef HeaderRaw
strdef HeaderExtras
strdef HeaderSomaArea
strdef HeaderShrinkageCorrection
strdef HeaderVersionNumber
strdef HeaderVersionDate
strdef HeaderScale
ScaleX = ScaleY = ScaleZ = 1
Shrinkagex = Shrinkagey = Shrinkagez = 1
SectionNum = 1

//------------------------------------------------------------------------------
// Section lists
//------------------------------------------------------------------------------
objref all, secInput, secOutput

//------------------------------------------------------------------------------
// synaptic points
//------------------------------------------------------------------------------
objref SynList
objref SynNo

objref synlist

objref synlist_
//------------------------------------------------------------------------------
// geometry
//------------------------------------------------------------------------------
tranlationX = translationY = translationZ = 0

//------------------------------------------------------------------------------
// procedures / functions interface
//------------------------------------------------------------------------------
// proc init()
// proc shape()
// proc translation()
// proc subsets()
// proc geom()
// proc biophys()
// proc geom_nseg()
// proc biophys_inhomo()
// obfunc makeSyn()
// obfunc makeNetCon

//
proc init() {
    
    SynList = new List()
    synlist = new List()
    SynNo = new Vector()
    
    synlist_ = new List()

    all = new SectionList()
    secInput = new SectionList()
    secOutput = new SectionList()
    
    shape($s1)
    subsets()
    geom()
    biophys()
    geom_nseg()
}

objref ID, Type, pID //pID -> parent ID
proc shape() { localobj strFunc, fobj strdef str, strName, strValue, strTemp
	fobj = new File()
	{fobj.ropen($s1)}
	strFunc = new StringFunctions()

	//--------------------------------------------------------------------------
	// load header information
	//--------------------------------------------------------------------------
	curLine = 0
	dataLineOverHead = 0
	while (!fobj.eof()) {
		fobj.gets(str)
		strFunc.head(str, " ", strName)
		strFunc.tail(strName, "#", strName)
		strFunc.tail(str, " ", strValue)
		strFunc.head(strValue, "\n", strValue)
		
		strFunc.head(str, "[^#]", strTemp)
		if (0 != strcmp("#", strTemp)) {
			dataLineOverHead = curLine
			break
		}
		
		if (!strcmp(strName, "ORIGINAL_SOURCE")) {
			HeaderOriginalSource = strValue
		}
		if (!strcmp(strName, "CREATURE")) {
			HeaderCreature = strValue
		}
		if (!strcmp(strName, "REGION")) {
			HeaderRegion = strValue
		}
		if (!strcmp(strName, "FIELD/LAYER")) {
			HeaderFieldLayer = strValue
		}
		if (!strcmp(strName, "TYPE")) {
			HeaderType = strValue
		}
		if (!strcmp(strName, "CONTRIBUTOR")) {
			HeaderContributor = strValue
		}
		if (!strcmp(strName, "REFERENCE")) {
			HeaderReference = strValue
		}
		if (!strcmp(strName, "RAW")) {
			HeaderRaw = strValue
		}
		if (!strcmp(strName, "EXTRAS")) {
			HeaderExtras = strValue
		}
		if (!strcmp(strName, "SOMA_AREA")) {
			HeaderSomaArea = strValue
		}
		if (!strcmp(strName, "SHRINKAGE_CORRECTION")) {
			HeaderShrinkageCorrection = strValue
			sscanf(HeaderShrinkageCorrection, "%f %f %f", &ShrinkageX, &ShrinkageY, &ShrinkageZ)
		}
		if (!strcmp(strName, "VERSION_NUMBER")) {
			HeaderVersionNumber = strValue
		}
		if (!strcmp(strName, "VERSION_DATE")) {
			HeaderVersionDate = strValue
		}
		if (!strcmp(strName, "SCALE")) {
			HeaderScale = strValue
			sscanf(HeaderScale, "%f %f %f", &ScaleX, &ScaleY, &ScaleZ)
			//printf("Scale (%f, %f, %f)\n", ScaleX, ScaleY, ScaleZ)
		}
		
		curLine += 1
	}
	
	SectionNum = 0
	while (!fobj.eof()) {
		SectionNum = fobj.scanvar()-1
		for i = 2, 7 {
			fobj.scanvar()
		}
	}
	
	{fobj.close()}
	
	{fobj.ropen($s1)}
	
	curLine = 0
	while (curLine < dataLineOverHead) {
	    fobj.gets(strTemp)
	    curLine += 1
	}
	
	//--------------------------------------------------------------------------------
	// create multi-compartment model
	//---------------------------------------------------------------------------------
	create Dend[SectionNum]
	ID = new Vector(SectionNum)
	Type = new Vector(SectionNum)
	pID = new Vector(SectionNum)

	// root node
	rootNo = fobj.scanvar() - 2
	type_root = fobj.scanvar()
	x_root = fobj.scanvar()
	y_root = fobj.scanvar()
	z_root = fobj.scanvar()
	d_root = fobj.scanvar()
	parentNo_root = fobj.scanvar() - 2
	
	while (1) { // while current node's connected to the root node
	    curNo = fobj.scanvar() - 2
	    type = fobj.scanvar()
	    x = fobj.scanvar()
	    y = fobj.scanvar()
	    z = fobj.scanvar()
	    d = fobj.scanvar()
	    parentNo = fobj.scanvar() - 2
	    ID.x[curNo]=curNo
	    Type.x[curNo]=type
	    pID.x[curNo]=parentNo
	    if (parentNo != rootNo){break}
	    Dend[curNo] {
		pt3dclear()
		pt3dadd(x_root * ScaleX, y_root * ScaleY, z_root * ScaleZ, d_root)
		pt3dadd(x * ScaleX, y * ScaleY, z * ScaleZ, d)
	    }
	}
	
	Dend[parentNo] {
	    parentX = x3d(1)
	    parentY = y3d(1)
	    parentZ = z3d(1)
	    parentDiam = diam3d(1)
	}
	Dend[curNo] {
	    pt3dclear()
	    pt3dadd(parentX, parentY, parentZ, parentDiam)
	    pt3dadd(x * ScaleX, y * ScaleY, z * ScaleZ, d)		
	}
	connect Dend[curNo](0), Dend[parentNo](1)
	
	while (!fobj.eof()) {
	    curNo = fobj.scanvar() - 2
	    type = 	fobj.scanvar()
	    x = fobj.scanvar()
	    y = fobj.scanvar()
	    z = fobj.scanvar()
	    d = fobj.scanvar()
	    parentNo = fobj.scanvar() - 2
	    ID.x[curNo]=curNo
	    Type.x[curNo]=type
	    pID.x[curNo]=parentNo
	    	    
	    Dend[parentNo] {
		parentX = x3d(1)
		parentY = y3d(1)
		parentZ = z3d(1)
		parentDiam = diam3d(1)
	    }
	    Dend[curNo] {
		pt3dclear()
		pt3dadd(parentX, parentY, parentZ, parentDiam)
		pt3dadd(x * ScaleX, y * ScaleY, z * ScaleZ, d)		
	    }
	    connect Dend[curNo](0), Dend[parentNo](1)
	}
	
	{fobj.close()}
	
	print "load : ", $s1, " ok"
} // shape()

proc translation() {
	translationX = $1
	translationY = $2
	translationZ = $3
	forall {
		x0 = x3d(0)
		y0 = y3d(0)
		z0 = z3d(0)
		diam0 = diam3d(0)
		x1 = x3d(1)
		y1 = y3d(1)
		z1 = z3d(1)
		diam1 = diam3d(1)
		pt3dclear()
		pt3dadd(translationX + x0, translationY + y0, translationZ + z0, diam0)
		pt3dadd(translationX + x1, translationY + y1, translationZ + z1, diam1)
	}
}

proc subsets() {
	Dend[0] all.wholetree()
}

proc geom() {
}

proc biophys() {
	// forall {
	// 	Ra = 100
	// 	cm = 1
	// }
}

external lambda_f
proc geom_nseg() {
	forall { nseg = int((L / (0.1 * lambda_f(100)) + .9) / 2) * 2 + 1 }
}

proc biophys_inhomo() {}

objref syn_
obfunc makeSyn() { local targetno localobj r
	// args : targetSectionNo, tau1, tau2, e
	for i = 0, SynNo.size-1 {
		if ($1 == SynNo.x[i]) return SynList.object(i)
	}
	if ($1 == -1) {
		r = new Random()
		r.uniform(0, SectionNum-1)
		targetno = int(r.repick())
	} else {
		targetno = $1
	}
	SynNo.append(targetno)
	Dend[targetno] syn_ = new Exp2Syn(0.5) SynList.append(syn_)
	syn_.tau1 = $2
	syn_.tau2 = $3
	syn_.e = $4
	return syn_
}
numSynapse = 0
objref SynNumberList
strdef SynF //Synapse List File name
proc ReadSynlist(){local i localobj f
    f = new File()
    SynNumberList = new Vector()
    SynF = $s1
    f.ropen(SynF)
    numSynapse = f.scanvar()
    for i=0, numSynapse-1{
	SynNumberList.append(f.scanvar())
    }
    //print numSynapse
    print SynNumberList
}

objref syn__

//fixed 2015.01.08
proc synapses(){local targetno
    targetno = $1
    Dend[targetno] syn__ = new ExpSid(0.2) synlist.append(syn__)
    //print syn__.e
    //syn__.tau = 2
    Dend[targetno] syn__ = new ExpSid(0.4) synlist.append(syn__)
    Dend[targetno] syn__ = new ExpSid(0.6) synlist.append(syn__)
    Dend[targetno] syn__ = new ExpSid(0.8) synlist.append(syn__)
    //syn__.tau = 5
    //syn__.e = -80
}

/*
proc synapse(){local targetno
    targetno = $1
    //    print "targetno ",targetno
    Dend[targetno] syn__ = new ExpSid(0.5) synlist.append(syn__)
    //    syn__.tau = 2
    //    syn__.e = -80
}
*/

objref Syn__
obfunc mkSyn(){local targetno
    targetno = $1
    
    Dend[targetno] Syn__ = new ExpSid(0.5) synlist_.append(Syn__) 
    Syn__.tau = 0.1
    Syn__.e   = $2
    return Syn__
}

public connect2target
obfunc connect2target(){local numSynapses,pre_comp localobj nc,target
    pre_comp = $1
    target = $o2
    //numSynapses = $3
    Dend[pre_comp] nc = new NetCon(&v(0.5),target)
    nc.threshold = 1
    //nc.weight = 10/numSynapses
    //print "nc.weight ",nc.weight
    return nc
}


public doNothing
proc doNothing(){
    print $s1
    print "do nothing"
}
strdef str1

//make synapses to connect RN and other cells
proc mkSynapses(){local i
    //$s1 = Synapse list data file name
    str1 = $s1
    ReadSynlist(str1)
    for i=0, SynNumberList.size()-1{
	synapses(SynNumberList.x[i])
    }
}

objref iclamp
proc mkIclamp(){local id, num,targetno
    targetno = $1
    //print "setSinglePulse", id
    Dend[targetno] iclamp = new IClamp(0.5)
    iclamp.dur = 20
    iclamp.amp = 1.5
    iclamp.del = 0.5
    //iclamp.del = 2*id*(num+1)
    //iclamp.del = 20*id*(num+1)
    //    psection()
}


obfunc makeNetCon() { local targetno localobj nc, r
	// args : sourceSectionNo, targetSyn, threshold, delay, weight
	if ($1 == -1) {
		r = new Random()
		r.uniform(0, SectionNum)
		targetno = int(r.repick())
	} else {
		targetno = $1
	}
	// NetCon(source, target, threshold, delay, weight)
	Dend[targetno] nc = new NetCon(&v(0.5), $o2, $3, $4, $5)
	return nc
}

endtemplate CellSwc