#pragma rtGlobals=1		// Use modern global access method.
//This procedure is an Igor procedure used to simulate the calcium clearance dynamics after 3-s depolarization. 
//Typing doall() in the command window in Igor Pro 4.0 will excute the whole procedure
Function DoAll()
	Make /o/N=500 store1=0, store2=0, store3=0, store4=0,store5=0,store6=0
	cc(0)
	cc(2)
	cc(3)
	cc(4)
	cc(5)
	cc(6)
	cc(7)
	Display/w=(10, 10,380,300)  store6, store5, store4, store3 vs TAxis;DelayUpdate 
	//store6=Na7.4,store5= Na8.8 or NaLa7.4
	//store4=Li7.4,store3= Li8.8 or LiLa7.4
	SetAxis bottom 0,100
	ShowInfo
	graph()
end function

//---------------------------------------------------GET Ca-BINDING RATIO OF BUFFER
Function GetBuffer(ci)
	variable ci
	variable B, cIndo = 0.188e-3 ,  kIndo = 0.0002e-3    // M units for indo conc and Ca
	B = 1+ 100+ (cIndo/kIndo)/(1+ci/kIndo)^2
	return B
end function

//---------------------------------------------------RATE CONSTANTS
Function Initialize Rates ()
//			all concentrations in M units
	variable BScale =1
	variable /g  kNCE=86.47e-6 *BScale
	variable /g  VMaxpmcaFast =21.02e-6*BScale
	variable /g VMaxpmcaSlow = VMaxpmcaFast/7.2   // Na8.8
	//variable /g VMaxpmcaSlow = VMaxpmcaFast/2.6  // NaLa7.4
	variable /g Kserca =270e-9, Kpmca=5e-7
	variable /g VMaxserca=0  //SERCA off
	//variable /g VMaxserca=40.92e-6*BScale  //SERCA on
	variable /g KmitoMCU =0* BScale, MitoBeg = 500e-9  // Mito off
	variable /g kRestLeak =2.2845e-4*BScale, kBigLeak = 349.13*kRestLeak//mM per sec
	variable /g NaoHS =130e-3, Erest = -70, Nai = 10e-3
end function	

//---------------------------------------------------Set PUMPS on or off
Function Adjust Rates (choice)
	Variable choice
// 0	control
// 1	Brief depolarization period
// 2 PMCA slowed (Na8.8 or NaLa7.4)
// 3	NCX inhibition (Li7.4)
// 4	PMCA slowed and NCX inhibition (Li8.8 or LiLa7.4)
// 5	PMCA slowed, NCX inhibition, Mito off
// 6 Just PMCA slowed, no depol
// 7 Mito off 

	variable /g VMaxPMCA, VMaxpmcaFast , VMaxpmcaSlow 
	variable /g kNCE
	variable /g KmitoMCU, MitoBeg, kMito
	variable /g kRestLeak, kBigLeak, kLeak
	variable /g  Nao, NaiSet,  NaoHS, Em, Erest
	
	//standard conditions
	kLeak = kRestLeak
	VMaxPMCA = VMaxPMCAFast
	Nao = NaoHS
	kMito = kMitoMCU
	Em = Erest
	NaiSet =10e-3
	variable /g dt=0.5, t
	Switch (choice)
		case 1:		//Kdepol, pH 7.4
			kLeak = kBigLeak* (0.60 + 1.8*exp(-(t-30)*2)) 
			Nao = 67.5e-3
			NaiSet = 0
			Em = -10
			dt=0.5
			break
		case 2:		// 2 PMCA slow (Na8.8 or NaLa7.4)
			VMaxPMCA = VMaxPMCASlow
			Kleak=kRestLeak
			kLeak = 1.6*kRestLeak
			break
		case 3:		// 3	NCE off (Li7.4)
			Nao = 0
			NaiSet = 0
			break 
		case 4:		// 4	PMCA slow and NCE off  (Li8.8 or NaLa7.4)
			VMaxPMCA = VMaxPMCASlow
			Kleak=kRestLeak
			kLeak = 1.6*kRestLeak
			Nao = 0
			NaiSet = 0			
			break
		case 5:		// 5	PMCA slow, NCE off, Mito off
			VMaxPMCA = VMaxPMCASlow
			Nao = 0
			NaiSet = 0			
			KMito = 0			
			break
		case 6:		//6 Just PMCA slow, no depol
			VMaxPMCA = VMaxPMCASlow
			kLeak =1.3*kRestLeak
			break
		case 7:		//7 Mitos off
			kMito=0
			break
							
	endswitch
end function

//--------------------------------------------------- CALCULATE FLUX from RATE LAWS
Function Flux()
	variable /g VMaxPMCA,kPMCA, kNCE, Kserca,Vmaxserca 
	variable /g NCE, PMCA, Mito, SERCA
	variable/g MitoBeg, kMito
	variable 	/g kLeak
	variable /g  Nao, Nai,Cai, Em
	variable/g Cao 
	NCE = - DB1(Cai,Cao,nai,nao,Em)
	PMCA = Vmaxpmca *(1/(1+(Kpmca/Cai)))
	SERCA =Vmaxserca*(1/(1+(Kserca/Cai)*(Kserca/Cai)))
	if (Cai < MitoBeg)
		Mito = 0
	else
		Mito = Kmito*(Cai-MitoBeg)
	endif
end function

/---------------------------------------------------INTEGRATE Cai TIME COURSE
Function cc(choice) //CellCalcium time course
	variable choice
//     start by declaring variables
	variable/g  VMaxPMCA, Kpmca
	variable/g Kmito, MitoBeg
	variable 	/g  KLeak
	variable /g kNCE
	variable /g  Nao,  CaiStart = 80e-9, Cai = CaiStart
	variable /g Cao = 2e-3
	variable dCaim
	variable /g NaiSet, dNai, Nai
	variable /g  t = -200, dt  = 0.5,  tEnd = 500
	variable i=t/dt
	variable/g NCE, Leak,PMCA,Mito, SERCA
	wave store1, store2, store3, store4,store5,store6    
	make /o /N=500 Camwave, TAxis, Naiwave
	store6 = store5
	store5 = store4
	store4 = store3
	store3 = store2
	store2 = store1
	store1 = Camwave
	InitializeRates()
	AdjustRates(0)
	
//    set initial values in wave
	TAxis=NaN
	TAxis[i]=t
	Camwave[i]= Cai

//      integrate Ca by Euler method until it reaches the End
	do
		if ((t>=30) && (t<33))//(i==3) 
			if (choice != 6)
				AdjustRates(1)
			End if
		Else if (t==33)//(i==5)
			Adjust Rates(choice)
			End if
		flux()
   		Leak = kLeak * Cao
		dCaim = (LEAK - PMCA - NCE - Mito-SERCA) * dt / GetBuffer(Cai)

//  integrate Nai as well		
		Cai += dCaim
		dNai = (3* NCE - 0.1*(Nai-NaiSet))*dt   //10 secs of equlibration time for Na
		Nai += dNai
		i += 1
		t  += dt
		Camwave[i] = Cai
		NaiWave[i]  = Nai
		TAxis[i] = t
	
	while (t<tEnd)
//	Scale stored values to nM concentration scale
	variable CScale = 1e9
	Camwave *= CScale
	
end function

//---------------------------------------------------[Ca} DEPENDENCE of rates
Function graph()			//DrawTransportRates()
//Draw a graph of each transport component
	//     start by declaring variables
	variable/g Cai =0
	variable dCaim=10e-9
	variable/g Cao = 2e-3  // outside in M
	variable da
	variable /g kNCE,VMaxPMCA, kLeak, kMito, kPMCA, MitoBeg
	variable /g  Nao, NaoHS
	Nao = NaoHS   
	variable i =1
	variable /g NCE, Leak,PMCA, Mito,SERCA
	InitializeRates()
	AdjustRates(0)
	Make /o/N=1000 CamStep
	CamStep=NaN
	make /o /N=200 LeakRate, NCERate, PMCARate, MitoRate, AllFlux, SERCARate
	LeakRate=NaN; NCERate=NaN; PMCARate=NaN; MitoRate=NaN;AllFlux=NaN
	variable CaiEnd = 2000e-9
	variable index =0
	//      step through values of Cai to make graph waves
	do
		leak = kleak*Cao
		Flux()
		Cai += dCaim
		LeakRate[index] = -Leak
		NCERate[index] = NCE
		PMCARate[index] = PMCA
		SERCARate[index] = SERCA

		MitoRate[index] = Mito
		Camstep[index] =Cai
		AllFlux[index] = NCE + PMCA + Mito - leak+SERCA
		index +=1
	
	while (Cai<CaiEnd)
	
	CamStep *= 1e9    //convert Calcium horizontal scale to nanomolar
	Display/w=(400, 10,750,400) LeakRate, NCERate, PMCARate, SERCARate, AllFlux  vs CamStep
	
end function

//---------------------------------------------------BERS models of NCX
function DB1(ci,co,ni,no,Em)
	//Don Bers empirical model for cardiac NCX current
	variable ci,ni,co,no,em
	variable kmnao=87.5e-3, kmnai=12.3e-3, kmcai=0.0036e-3, kmcao=1.3e-3, ksat=0.27, nu=0.35
	variable inaca, vmax=0.8, rt=25
	variable /g kNCE
	inaca  = kNCE*vmax*(ni^3*co*exp(nu*Em/rt) - no^3*ci*exp((nu-1)*Em/rt))
	inaca /= (kmcao*ni^3+kmnao^3*ci+kmnai^3*co*(1+(ci/kmcai))+kmcai*no^3*(1+(ni/Kmnai)^3)+ni^3*co+no^3*ci)
	inaca /= (1+ksat*exp((nu-1)*Em/rt))
	return inaca
end function