CA1PC.hoc

// General CA1 pyramidal cell model
//   - with minimal set of ion channels and their distributions
// Morphology of cell n5038804 as used in Migliore et al J Neurophys 94:4145-4155, 2005.
// Based on Migliore et al 2005 - passive properties, Na, Kdr, Ka, Ih 
// and Poirazzi et al Neuron 37:977-987, 2003 - IM, calcium currents (L, T, R), mAHP, sAHP
// Also Holbro et al PNAS 107:15975-15980, 2010 (fast Ca R)

// Geometry files must define the section lists:
// basal_list, apical_list, trunk_list, oblique_list, SR_list, SLM_list
// soma and axon sections must be called "soma" and "axon"
// SLMbr_list is a special list containing the sections in a single SLM branch
// (used for focal stimulation)

// Last update: BPG 2-5-14

begintemplate PyramidalCell
public is_art, totnseg, totarea, Vrest
public gna, NaMULT, arSOMA, arDIST, arMAX, gka, ghd
public init, topol, basic_shape, subsets, geom, biophys, set_dendrite
public pre_list, connect2target
public cell_init, plotvmx, newsyn, showsyn, showlayersyn, plotcamx, plotcamxd, plotcamxt

public soma, soma_list, axon_list, dendrite_list
public basal_list, apical_list, trunk_list, oblique_list
public SR_list, SRbrp_list, SRbrd_list, SRbrt_list, SRbrt2_list, SRbrtc_list
public SRprox_list, SRdist_list
public SLM_list, SLMbr_list, SLMbrt_list
public shead, sneck, spine_list, sneck_list

objref basal_list, apical_list, trunk_list, oblique_list
objref SR_list, SRbrp_list, SRbrd_list, SRbrt_list, SRbrt2_list, SRbrtc_list
objref SRprox_list, SRdist_list
objref SLM_list, SLMbr_list, SLMbrt_list
objref soma_list, axon_list, dendrite_list, pre_list, spine_list, sneck_list


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


xopen("Cells/Mig_sec5038804.hoc")             // Migliore section file


proc init() {

soma_list = new SectionList()
axon_list = new SectionList()
dendrite_list = new SectionList()
basal_list = new SectionList()
apical_list = new SectionList()
trunk_list = new SectionList()
oblique_list = new SectionList()
SR_list = new SectionList()
SRbrp_list = new SectionList()
SRbrd_list = new SectionList()
SRbrt_list = new SectionList()
SRbrt2_list = new SectionList()
SRbrtc_list = new SectionList()
SRprox_list = new SectionList()
SRdist_list = new SectionList()
SLM_list = new SectionList()
SLMbr_list = new SectionList()
SLMbrt_list = new SectionList()
spine_list = new SectionList()
sneck_list = new SectionList()

xopen("Cells/Mig_geo5038804.hoc")             // Migliore geometry file
fix_seg = 0

// Spine dimesions (if spines used)
// approximate Grunditz et al model
//sneck_diam = 0.04 /* um (high resistance) */
sneck_diam = 0.125 /* um (low resistance) */
sneck_len = 1.0 /* um */
shead_diam = 0.5 /* um */
shead_len = 0.5 /* um */

// Passive properties from Migliore 2005
Vrest = -65
celsius = 35.0  

Rm = 28000	// Migliore
//Rm = 200000	// Poirazi soma (cf Kali 300000)
RmDend = Rm
RmSoma = Rm
RmAx = Rm

Cm    = 1
CmSoma= Cm
CmAx  = Cm
CmDend = Cm

//RaAll= 50	// low Ra
//RaSoma=50  
//RaAx = 50
//RaSpine = 50
RaAll= 150	// Migliore 2005
RaSoma=150  
RaAx = 150
RaSpine = 150

gna =  0.025	// Migliore 2005
AXONM = 5
NaMULT = 0.015/0.025	// scale Na in dendrites (=1 for Migliore)
arSOMA = 1.0	// slow inactivation at soma (1=none; 0=max) (BPG)
arMAX = 0.5	// slow inact at distance (BPG)
arDIST = 350	// distance at which inactivation saturates (BPG)

gkdr = 0.01

gka = 0.03	// standard high KA
//gka = 0.01	// suitably reduced: low KA
KMULT =  gka
KMULTP = gka

ghd=0.00005
//ghd=0.0001	// doubled Ih

soma_km = 0.06	// not used (BPG)

soma_caR = 0.03 
soma_sAHP = 0.001
soma_mAHP = 0.001

caR_init = 0.03 
sAHP_init = 0.001
mAHP_init = 0.001
caR_spine = 0.03 
sAHP_spine = 0.001
mAHP_spine = 0.001
caR_vact = -30	// make half-activation point more hyperpolarised
caR_tinact = 20	// make inactivation faster


forsec axon_list {insert pas e_pas=Vrest g_pas = 1/RmAx Ra=RaAx cm=CmAx}
forsec soma_list {insert pas e_pas=Vrest g_pas = 1/RmSoma Ra=RaSoma cm=CmSoma}
forsec dendrite_list {insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}

access soma
freq=50
if (!fix_seg) {geom_nseg()}	// adjust segment number in each section
totnseg=0
forall {totnseg=totnseg+nseg}
totarea=0
forall {for (x) totarea=totarea+area(x)}
access soma
distance()

forsec axon_list {   
                insert nax gbar_nax=gna*AXONM
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
}

forsec soma_list {   
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
                insert na3 gbar_na3=gna ar_na3=arSOMA
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
       		//insert km gbar_km=soma_km
            	insert carF gcabar_carF = soma_caR vha_carF=caR_vact ti_carF=caR_tinact           
            	//insert kca gbar_kca = soma_sAHP  // sAHP       
            	insert kmAHP gkbar_kmAHP = soma_mAHP  // mAHP      
       		insert cad    // calcium pump/buffering mechanism
}

forsec dendrite_list {
		insert hd ghdbar_hd=ghd
                insert na3 gbar_na3=gna*NaMULT
                insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0
       		//insert km gbar_km=soma_km
            	insert carF gcabar_carF = caR_init vha_carF=caR_vact ti_carF=caR_tinact            
            	//insert kca gbar_kca = sAHP_init  // sAHP       
            	insert kmAHP gkbar_kmAHP = mAHP_init  // mAHP      
       		insert cad    // calcium pump/buffering mechanism
       		insert ds
       		insert dca

		for (x) if (x>0 && x<1) { xdist = distance(x)
		
			// Na slow inactivation saturates at 350um (BPG)
			if (xdist < arDIST) {
                	  ar_na3(x) = arSOMA-(arSOMA-arMAX)*(xdist/arDIST)
                	} else {
                	  ar_na3(x) = arMAX
                	}

			// Ih saturates at 350um (BPG)
			if (xdist < 350) {
                	  ghdbar_hd(x) = ghd*(1+3*xdist/100)
                	} else {
                	  ghdbar_hd(x) = ghd*(1+3*350/100)
                	}
                	
                	// KA saturates at 350um (BPG)
                	if (xdist > 100 && xdist < 350) {
				vhalfl_hd=-81
                        	gkabar_kad(x) = KMULT*(1+xdist/100)
                        } else if (xdist >= 350) {
				vhalfl_hd=-81
                        	gkabar_kad(x) = KMULT*(1+350/100)
                	} else {
				vhalfl_hd=-73
                        	gkabar_kap(x) = KMULTP*(1+xdist/100)
               		}
		}
}

cell_init()

pre_list = new List()
//synapses()


} // end init


proc set_dendrite() {
	KMULT =  gka
	KMULTP = gka
	access soma
	distance()
	forsec dendrite_list {
		gbar_na3=gna*NaMULT		
		for (x) if (x>0 && x<1) { xdist = distance(x)
		
			// Na slow inactivation saturates at 350um (BPG)
			if (xdist < arDIST) {
                	  ar_na3(x) = arSOMA-(arSOMA-arMAX)*(xdist/arDIST)
                	} else {
                	  ar_na3(x) = arMAX
                	}

			// Ih saturates at 350um (BPG)
			if (xdist < 350) {
                	  ghdbar_hd(x) = ghd*(1+3*xdist/100)
                	} else {
                	  ghdbar_hd(x) = ghd*(1+3*350/100)
                	}
                	
                	// KA saturates at 350um (BPG)
                	if (xdist > 100 && xdist < 350) {
				vhalfl_hd=-81
                        	gkabar_kad(x) = KMULT*(1+xdist/100)
                        } else if (xdist >= 350) {
				vhalfl_hd=-81
                        	gkabar_kad(x) = KMULT*(1+350/100)
                	} else {
				vhalfl_hd=-73
                        	gkabar_kap(x) = KMULTP*(1+xdist/100)
               		}
		}
	}
	cell_init()
}


proc cell_init() {
	t=0
        forall {
        v=Vrest
        if (ismembrane("nax") || ismembrane("na3")) {ena=55}
        if (ismembrane("kdr") || ismembrane("kap") || ismembrane("kad")) {ek=-90}
        if (ismembrane("hd") ) {ehd_hd=-30}
	}
	finitialize(Vrest)
        fcurrent()

        forall {
	for (x) {
	if (ismembrane("na3")||ismembrane("nax")){e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)}
	if (ismembrane("hd")) {e_pas(x)=e_pas(x)+i_hd(x)/g_pas(x)}
		}
	}
}


objref distrx, distry, c
proc plotvmx() {	// plot max voltage with distance
	c = new Graph()
	c.size(0,1000,0,100)
	c.xaxis(1)
	c.exec_menu("10% Zoom out")
	c.color(1)
	c.label(0.4,0.8," peak AP")

	distrx=new Vector()
	distry=new Vector()
	forsec $o1 {
		for (x) if (x>0 && x<1) {
			if (diam>=0.) {
			distrx.append(distance(x)) 
			distry.append(vmax_ds(x)-Vrest)
			}
		}
	}
	distry.mark(c,distrx,"O",3,3,2)
//	distry.mark(c,distrx,"t",5,1,1) 
	c.flush()
	doNotify()
}


objref cax, cay, cca
proc plotcamx() {local i	// plot max ca across spine heads
	cca = new Graph()
	cca.size(0,$2,0,$3)
	cca.xaxis(1)
	cca.exec_menu("10% Zoom out")
	cca.color(1)
	cca.label(0.4,0.8," peak Ca")

	cax=new Vector()
	cay=new Vector()
	i = 0
	forsec $o1 {
	  i = i+1
	  cax.append(i) 
	  cay.append(camax_dca(0.5))
	}
	cay.mark(cca,cax,"O",3,3,2)
	cca.flush()
	doNotify()
	print cay.mean(), cay.max(), cay.min()
}

objref tcax, tcay, ccat
proc plotcamxt() {local i	// plot max ca against time
	ccat = new Graph()
	ccat.size(0,$2,0,$3)
	ccat.xaxis(1)
	ccat.exec_menu("10% Zoom out")
	ccat.color(1)
	ccat.label(0.4,0.8," peak Ca")

	tcax=new Vector()
	tcay=new Vector()
	i = 0
	forsec $o1 {
	  tcax.append(tmax_dca(0.5)) 
	  tcay.append(camax_dca(0.5))
	}
	tcay.mark(ccat,tcax,"O",3,3,2)
	ccat.flush()
	doNotify()
	print tcay.mean(), tcay.max(), tcay.min()
}

objref dcax, dcay, ccad
proc plotcamxd() {	// plot max calcium with distance
	ccad = new Graph()
	ccad.size(0,1000,0,1)
	ccad.xaxis(1)
	ccad.exec_menu("10% Zoom out")
	ccad.color(1)
	ccad.label(0.4,0.8," peak Ca")

	dcax=new Vector()
	dcay=new Vector()
	forsec $o1 {
		for (x) if (x>0 && x<1) {
			if (diam>=0.) {
			dcax.append(distance(x)) 
			dcay.append(camax_dca(x))
			}
		}
	}
	dcay.mark(ccad,dcax,"O",3,3,2)
//	distry.mark(ccad,dcax,"t",5,1,1) 
	ccad.flush()
	doNotify()
	print dcay.mean(), dcay.max(), dcay.min()
}


obfunc connect2target() { localobj nc //$o1 target point process, optional $o2 returned NetCon
  soma nc = new NetCon(&v(1), $o1)
  nc.threshold = -10
  if (numarg() == 2) { $o2 = nc } // for backward compatibility
  return nc
}


// Code to add individual synapses in random locations (BPG 6-2-10)
create shead[1], sneck[1]
spi = 0		// current spine index

proc newsyn() { local i, j localobj syn, rl, syn_list, snr, shr 
// $1 type $2 number $o3 target section list, $o4 uniform random no., $5 flag spines, $6 total spines
  access soma
  rl = new RandomLocation($o3, $o4)
  if ($5 == 1 && spi == 0) {
    create sneck[$6], shead[$6]
  }
  for i=0, $2-1 {
    syn_list = new List()
    if ($1 == 1) {	// AMPA
      soma syn = new Exp2Syn(0.5) pre_list.append(syn)
      syn.tau1 = 0.5
      syn.tau2 = 3
      syn.e = 0
      syn_list.append(syn)
    } else if ($1 == 2) {	// AMPA/NMDA
      soma syn = new Exp2Syn(0.5) pre_list.append(syn)
      syn.tau1 = 0.5	// Spruston JPhysiol 1995
      syn.tau2 = 3	// Spruston JPhysiol 1995
      syn.e = 0
      syn_list.append(syn)
      soma syn = new NMDAca(0.5) pre_list.append(syn)
      syn.fCa = 0.1	// fraction of Ca current (Bloodgood & Sabatini)
      syn.tcon = 3	
      syn.tcoff = 100
      syn.mgconc = 1	// (mM) standard Mg conc
      syn.gamma = 0.08	// Larkum Science 2009 (sharpens voltage curve)
      syn_list.append(syn)
    } else if ($1 == 3) {	// GABAA
      soma syn = new Exp2Syn(0.5) pre_list.append(syn)
      syn.tau1 = 1
      syn.tau2 = 8
      syn.e = -75
      syn_list.append(syn)
    } else if ($1 == 4) {	// GABAB
      soma syn = new Exp2Syn(0.5) pre_list.append(syn)
      syn.tau1 = 35
      syn.tau2 = 100
      syn.e = -75
      syn_list.append(syn)
    }
    if ($5 == 1) {  // spines
      sneck[spi].L = sneck_len
      sneck[spi].diam = sneck_diam
      shead[spi].L = shead_len
      shead[spi].diam = shead_diam
      sneck[spi] {insert pas e_pas=Vrest g_pas=1/RmDend Ra=RaAll cm=CmDend
        insert cad taur_cad=14 insert ds insert dca}
      shead[spi] {insert pas e_pas=Vrest g_pas=1/RmDend Ra=RaAll cm=CmDend
        insert carF gcabar_carF=caR_spine vha_carF=caR_vact ti_carF=caR_tinact
        insert kmAHP gkbar_mAHP = mAHP_spine  // mAHP      
        insert cad taur_cad=14 depth_cad=shead_diam/2
        insert ds insert dca 
        }
      sneck[spi] snr = new SectionRef()
      sneck[spi] sneck_list.append()
      shead[spi] shr = new SectionRef()
      shead[spi] spine_list.append()
      access soma
      rl.locsp(syn_list, snr, shr)
      spi = spi + 1
    } else {	// no spines
      rl.loc(syn_list)
    }
  }
}

objref cell_shape

proc showsyn() { local i
  cell_shape = new Shape()
  for i=0,pre_list.count()-1 {
    cell_shape.point_mark(pre_list.o(i), 3)
  }
}

proc showlayersyn() { local i
  if ($4 == 1) {
    cell_shape = new Shape()
  }
  for i=$1,$2 {
    cell_shape.point_mark(pre_list.o(i), $3)
  }
}


//**********************************************************************

objref syn_
proc synapses() {
  	/* E0 */   	soma syn_ = new Exp2Syn(0.5)  pre_list.append(syn_)
    	syn_.tau1 = 0.5
    	syn_.tau2 = 3
}


func is_art() { return 0 }

endtemplate PyramidalCell