: Intracellular calcium model 
: the calcium flux is in mM/ms

NEURON {
SUFFIX caodmy
USEION c WRITE ci VALENCE 0
USEION myca READ imyca WRITE mycai VALENCE 2
USEION ca READ ica
:RANGE ton1,toff1,ton2,toff2,ton3,toff3,ton4,toff4,ton5,toff5,conc1,conc2,conc3,ao
: RANGE ton, toff, ton1, toff1, conc1, conc2, ao, fca1, fca2
RANGE del, dur, ao, fca1, fca2
RANGE vrvgcc,gbar
GLOBAL jvgcc, RateLa2, caintra
}

UNITS {
	(mV) = (millivolt)
	(mA) = (milliamp)
	FARADAY = (faraday) (coulombs)
	PI = (pi) (1)
	(molar) = (1/liter)
	(mM)=(millimolar)
}



PARAMETER {

    :convfact=1 (cm2-mM/ms-mA)
    :convfact=3.79 (cm2-mM/ms-mA)
    :convfact2=1e-6 (mM/ms-mV)

	vol=31.16e-12 (cm3)
	surf=238.18e-8 (cm2)
	
	fca1=1e-3
	fca2=1e-3
	Fc=96485 (coul)

	: parameters from Ususyama et al 2012
	ct0=1.0000000000000002e-3  (mM)
	ct1=-1.1102230246251565e-19 (mM)
	ct2=1e-3 (mM)
	ct3=10.000000000000002e-3 (mM)
	ct4=1.0000000000000002e-3 (mM)
    KCaMCaK=4.06 (1/mM-ms)
	KCaMCa=0.64e-3 (1/ms)
	Ks=0.5  (1/mM-ms)
	Kd=45.44e-3 (1/ms)
	Ks1=10.54 (1/mM-ms)
	Kd1=15.15e-3 (1/ms)
    k1=50.76 (1/mM-ms)
    k2=0.12e-3 (1/ms)
    k11=1.35 (1/mM-ms)
	k21=0.1e-3 (1/ms)
	k12=0.16 (1/mM-ms)
    k22=0.01e-3 (1/ms)
	k13=1.01 (1/mM-ms)
	k23=0.73e-3 (1/ms)
	k14=23.16 (1/mM-ms)
	k24=0.43e-3 (1/ms)
	k15=0.04 (1/mM-ms)
	k25=0.78e-3 (1/ms)
	k16=1.55 (1/mM-ms)
	k26=0.01e-3 (1/ms)
	v1=37.3e-6 (mM/ms)
	k17=4.98e-3 (1/ms)
	k18=0.139 (1/mM-ms)
	k27=0.02e-3 (1/ms)
	k19=0.01e-3 (1/ms)
	k110=0.36 (1/mM-ms)
	k28=0.01e-3 (1/ms)
	k111=0.012e-3 (1/ms)
	k112=1.1 (1/mM-ms)
	k29=0.01e-3 (1/ms)
	k113=25.78e-3 (1/ms)
	k114=1.73 (1/mM-ms)
	k210=7e-3 (1/ms)
	k115=0.01e-3 (1/ms)
	k116=19.69 (1/mM-ms)
	k211=0.01e-3 (1/ms)
	k117=36.41e-3 (1/ms)
	:EfuCNG=1 
	:EfuVGCC=1
	:VRuVGCC=60
	:V50uVGC=-0.4878
	:KuVGCC=9.6469
	EfuCax=67.51e-3 (1/ms)
	KuCaX=19.04e-3 (mM)
	nuCaX=0.1
	k118=30.48 (1/mM-ms)
	k212=500.61e-3 (1/ms)
	k119=6.48 (1/mM-ms)
	k213=1.7e-3 (1/ms)
	k120=0.1 (1/mM-ms)
	k214=0.04e-3 (1/ms)
	Rutotal=1e-3 (mM)
    Gutotal=1e-3 (mM)
	ton (ms)
	toff (ms)
	ton1 (ms)
	toff1 (ms)
	:ton2 (ms)
	:toff2 (ms)
	:ton3 (ms)
	:toff3 (ms)
	:ton4 (ms)
	:toff4 (ms)
	:ton5 (ms)
	:toff5 (ms)
    conc1 (mM)   
    conc2 (mM)
	del=80000 (ms)
	dur=10000 (ms)
	: freq=1 (Hz)
    :conc3 (mM)
    :conc4 (mM)  
    :vrvgcc = 30 (mV)
    vrvgcc = 60 (mV)
	kvgcc=9.6469 (mV)
    v50vgcc=-0.4878 (mV)
	:kvgcc=4 (mV)
    :v50vgcc=-10 (mV)
}


ASSIGNED {
		v (mV)
		ao (mM)
		:icacng (mA/cm2) 
		icavgcc (mA/cm2) 
		ica (mA/cm2) 
		imyca (mA/cm2) 
		GGCYGC (mM)
		GCAPbCa (mM)
		GCAPaCa (mM)
		CaM (mM)
		PDEcGMP (mM)
		RateLa1 (mM/ms)
		RateLa2 (mM/ms)
		RateLaw (mM/ms)
		RateLa (mM/ms)
		MassAct (mM/ms)
		MassAc (mM/ms)
		MassAc1 (mM/ms)
		MassAc2 (mM/ms)
		MassAc3 (mM/ms)
		MassAc4 (mM/ms)
		MassAc5 (mM/ms)
		Constan (mM/ms)
		MassAc6 (mM/ms)
		MassAc7 (mM/ms)
		MassAc8 (mM/ms)
		MassA (mM/ms)
		MassA1 (mM/ms)
		MassA2 (mM/ms)
		MassA3 (mM/ms)
		MassA4 (mM/ms)
		MassA5 (mM/ms)
		MassA6 (mM/ms)
		MassA7 (mM/ms)
		RateLa3 (mM/ms)
		MassA8 (mM/ms)
		MassA9 (mM/ms)
		MassA10 (mM/ms)
		MassA11 (mM/ms)
		MassA12 (mM/ms)
		MassA13 (mM/ms)
		MassA14 (mM/ms)
		jvgcc (mM/ms)
		caintra (mM)
		}


STATE { 
	mycai (mM)
	ci (mM)
	GTP	(mM)
	GCY	(mM)
	Guactiv (mM) 
	PDE (mM)
	GCYGCAP	(mM)
	GGCY	(mM)
	PDEuact (mM)
	GCAPb (mM)
	CaMCa (mM)
    CaMCa3 (mM)
	GCYGCA	(mM)
	GCAPa	(mM)
	CaMCa2	(mM)
	Ruactiv	(mM)
	GGCYGCA (mM)
    CaMCa4	(mM)
	GCYGCA1 (mM)
	GCYGTP  (mM)
	GCYGCA2	(mM)
	PDEuac	(mM)
	}

INITIAL{
 	mycai=0.06e-3
	ci=0.54e-3
	GTP=7.48e-3
	GCY=7e-7
	Guactiv=0
	PDE=0.87e-3
	GCYGCAP=7.8e-6	
	GGCY=0
	PDEuact=0.012e-3
	GCAPb=0.86e-3
	CaMCa=2.37e-3
    CaMCa3=0.344e-3
	GCYGCA=7.7e-6
	GCAPa=0.0270e-3
	CaMCa2=1.1e-3
	Ruactiv=1.1e-8
	GGCYGCA=4.014e-11
    CaMCa4=0.13e-3
	GCYGCA1=2.5e-6
	GCYGTP=0.023e-3
	GCYGCA2=0.958e-3
	PDEuac=2e-6
	ao = conc1
}


BREAKPOINT {
       SOLVE states METHOD cnexp
	   jvgcc=RateLa2
	   caintra=mycai
	   }


DERIVATIVE states {

: Odorant 
	if (t> del && t<=del+dur){
	
	ao = conc1+((t-del)*((conc2-conc1)/((del+dur)-del)))

	} else{
	 ao=conc1
	 
	}

RateLaw=Ks*ao*(Rutotal-Ruactiv)-Kd*Ruactiv
RateLa=Ks1*Ruactiv*(Gutotal-Guactiv)-Kd1*Guactiv
MassAct=k1*GCY*GCAPa-k2*GCYGCAP
MassAc=k11*GCY*GCAPb-k21*GCYGCA
MassAc1=k12*GCY*Guactiv-k22*GGCY
MassAc2=k13*GCYGCAP*Guactiv-k23*GGCYGCA
MassAc3=k14*GCYGCA*Guactiv-k24*GGCYGC
MassAc4=k15*GGCY*GCAPa-k25*GGCYGCA
MassAc5=k16*GGCY*GCAPb-k26*GGCYGC
Constan=v1
MassAc6=k17*GTP
MassAc7=k18*GCY*GTP-k27*GCYGTP
MassAc8=k19*GCYGTP
MassA=k110*GCYGCAP*GTP-k28*GCYGCA2
MassA1=k111*GCYGCA2
MassA2=k112*GCYGCA*GTP-k29*GCYGCA1
MassA3=k113*GCYGCA1
MassA4=k114*PDE*ci-k210*PDEcGMP
MassA5=k115*PDEcGMP
MassA6=k116*PDEuact*ci-k211*PDEuac
MassA7=k117*PDEuac
MassA8=KCaMCaK*CaM*mycai-KCaMCa*CaMCa
MassA9=KCaMCaK*CaMCa*mycai-KCaMCa*CaMCa2
MassA10=KCaMCaK*CaMCa2*mycai-KCaMCa*CaMCa3
MassA11=KCaMCaK*CaMCa3*mycai-KCaMCa*CaMCa4
MassA12=k118*PDE*CaMCa4-k212*PDEuact
MassA13=k119*GCAPa*mycai-k213*GCAPaCa
MassA14=k120*GCAPb*mycai-k214*GCAPbCa


:RateLa1=convfact*imyca
:RateLa2=convfact2*(v-vrvgcc)/(1+exp((v50vgcc-v)/kvgcc))

RateLa1=fca1*(imyca*surf)/(2*vol*Fc)
:RateLa1=fca1*(imyca*surf*1e-3)/(2*vol*Fc)
RateLa2=fca2*(ica*surf)/(2*vol*Fc)
RateLa3=EfuCax*mycai/(1+(CaMCa4/KuCaX)^nuCaX)


: metabolite 'GuGCYuGCAPb' reactions
GGCYGC=ct0-1*GCY-GCYGCAP-1*GGCY-1*GCYGCA-GGCYGCA-GCYGCA1-1*GCYGTP-1*GCYGCA2
: metabolite 'GCAPbuCa' reactions
GCAPbCa=ct1+1*GCY+GCYGCAP+1*GGCY-1*GCAPb+GGCYGCA+1*GCYGTP+1*GCYGCA2
: metabolite 'GCAPauCa' reactions
GCAPaCa=ct2-GCYGCAP-GCAPa-GGCYGCA-GCYGCA2
: metabolite 'CaM' reactions
CaM=ct3-PDEuact-CaMCa-CaMCa3-CaMCa2-CaMCa4-PDEuac
: metabolite 'PDEucGMP' reactions
PDEcGMP=ct4-PDE-PDEuact-PDEuac


:Equations
mycai'=(-RateLa1-RateLa2-RateLa3-MassA8-MassA9-MassA10-MassA11-MassA13-MassA14)
:mycai'=(-RateLa1-RateLa3-MassA8-MassA9-MassA10-MassA11-MassA13-MassA14)
:mycai'=(-RateLa3-MassA8-MassA9-MassA10-MassA11-MassA13-MassA14)
ci'=(MassAc8+MassA1+MassA3-MassA4-MassA6)
GTP'=(Constan-MassAc6-MassAc7-MassA-MassA2)
GCY'=(-MassAct-MassAc-MassAc1-MassAc7+MassAc8)
Guactiv'=(RateLa-MassAc1-MassAc2-MassAc3)
PDE'=(-MassA4+MassA5-MassA12)
GCYGCAP'=(MassAct-MassAc2-MassA+MassA1)
GGCY'=(MassAc1-MassAc4-MassAc5)
PDEuact'=(-MassA6+MassA7+MassA12)
GCAPb'=(-MassAc-MassAc5-MassA14)
CaMCa'=(MassA8-MassA9)
CaMCa3'=(MassA10-MassA11)
GCYGCA'=(MassAc-MassAc3-MassA2+MassA3)
GCAPa'=(-MassAct-MassAc4-MassA13)
CaMCa2'=(MassA9-MassA10)
Ruactiv'=(RateLaw)
GGCYGCA'=(MassAc2+MassAc4)
CaMCa4'=(MassA11-MassA12)
GCYGCA1'=(MassA2-MassA3)
GCYGTP'=(MassAc7-MassAc8)
GCYGCA2'=(MassA-MassA1)
PDEuac'=(MassA6-MassA7)
}