/*
 * Calcium and Glycolysis Mediate Multiple Bursting Modes in Pancreatic
 * Islets
 * 
 * Model Status
 * 
 * This model has been rebuilt and coded by translating the authors
 * original XPPAUT.ode file, which can be found at http://www.math.fsu.edu/~bertram/software/islet/BJ_04a.ode.
 * This file runs in PCEnv and COR and is able to produce the expected
 * output. This model has been parameterised for the 'compound
 * bursting' model. Some equations have inconsistent magnitudes
 * due to the fact that time is defined in milliseconds, yet flux
 * is defined in micromolar per second.
 * 
 * Model Structure
 * 
 * ABSTRACT: Pancreatic islets of Langerhans produce bursts of
 * electrical activity when exposed to stimulatory glucose levels.
 * These bursts often have a regular repeating pattern, with a
 * period of 10-60 s. In some cases, however, the bursts are episodic,
 * clustered into bursts of bursts, which we call compound bursting.
 * Consistent with this are recordings of free Ca2+ concentration,
 * oxygen consumption, mitochondrial membrane potential, and intraislet
 * glucose levels that exhibit very slow oscillations, with faster
 * oscillations superimposed. We describe a new mathematical model
 * of the pancreatic beta-cell that can account for these multimodal
 * patterns. The model includes the feedback of cytosolic Ca2+
 * onto ion channels that can account for bursting, and a metabolic
 * subsystem that is capable of producing slow oscillations driven
 * by oscillations in glycolysis. This slow rhythm is responsible
 * for the slow mode of compound bursting in the model. We also
 * show that it is possible for glycolytic oscillations alone to
 * drive a very slow form of bursting, which we call "glycolytic
 * bursting." Finally, the model predicts that there is bistability
 * between stationary and oscillatory glycolysis for a range of
 * parameter values. We provide experimental support for this model
 * prediction. Overall, the model can account for a diversity of
 * islet behaviors described in the literature over the past 20
 * years.
 * 
 * Calcium and Glycolysis Mediate Multiple Bursting Modes in Pancreatic
 * Islets, Richard Bertram, Leslie Satin, Min Zhang, Paul Smolen,
 * and Arthur Sherman, 2004, Biophysical Journal, 87, 3074-3087.
 * PubMed ID: 15347584
 * 
 * cell diagram
 * 
 * [[Image file: bertram_2004.png]]
 * 
 * A schematic diagram of the ionic currents and fluxes across
 * the ER and the cell surface membranes which are described by
 * the mathematical model.
 */

import nsrunit;
unit conversion on;
// unit millisecond predefined
// unit millivolt predefined
// unit micromolar predefined
unit nanoS=1E-9 kilogram^(-1)*meter^(-2)*second^3*ampere^2;
unit picoS=1E-12 kilogram^(-1)*meter^(-2)*second^3*ampere^2;
unit femtoF=1E-15 kilogram^(-1)*meter^(-2)*second^4*ampere^2;
unit femtoA=1E-15 ampere^1;
unit per_second=1 second^(-1);
unit first_order_rate_constant=1E3 second^(-1);
unit micromolar_per_femtoA_millisecond=1E15 meter^(-3)*second^(-1)*ampere^(-1)*mole^1;
unit flux=1E-3 meter^(-3)*second^(-1)*mole^1;

math main {
	realDomain time millisecond;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) millivolt;
	when(time=time.min) V=-60;
	real IK(time) femtoA;
	real ICa(time) femtoA;
	real IKCa(time) femtoA;
	real IKATP(time) femtoA;
	real Cm femtoF;
	Cm=5300;
	real gK picoS;
	gK=2700;
	real VK millivolt;
	VK=-75;
	real n(time) dimensionless;
	when(time=time.min) n=0;
	real gKCa picoS;
	gKCa=600;
	real kd micromolar;
	kd=0.5;
	real c(time) micromolar;
	when(time=time.min) c=0.25;
	real gCa picoS;
	gCa=1000;
	real minf(time) dimensionless;
	real VCa millivolt;
	VCa=25;
	real taun millisecond;
	taun=20;
	real ninf(time) dimensionless;
	real cer(time) micromolar;
	when(time=time.min) cer=185;
	real fcyt dimensionless;
	fcyt=0.01;
	real Jmem(time) flux;
	real Jer(time) flux;
	real fer dimensionless;
	fer=0.01;
	real sigmaV dimensionless;
	sigmaV=31;
	real pleak first_order_rate_constant;
	pleak=0.0002;
	real Kserca first_order_rate_constant;
	Kserca=0.4;
	real lambdaer dimensionless;
	lambdaer=1;
	real epser dimensionless;
	epser=1;
	real alpha micromolar_per_femtoA_millisecond;
	alpha=0.00000450;
	real kpmca first_order_rate_constant;
	kpmca=0.2;
	real was_Jserca(time) flux;
	real Jleak(time) flux;
	real rgpdh(time) flux;
	real Rgk per_second;
	Rgk=0.2;
	real atot micromolar;
	atot=3000;
	real pfkbas dimensionless;
	pfkbas=0.06;
	real f6p(time) micromolar;
	real lambda dimensionless;
	lambda=0.005;
	real pfk(time) micromolar;
	real g6p(time) micromolar;
	when(time=time.min) g6p=200;
	real fbp(time) micromolar;
	when(time=time.min) fbp=40;
	real bottom1 dimensionless;
	bottom1=1;
	real topa1 dimensionless;
	topa1=0;
	real k1 micromolar;
	k1=30;
	real k2 micromolar;
	k2=1;
	real k3 micromolar;
	k3=50000;
	real k4 micromolar;
	k4=1000;
	real cat dimensionless;
	cat=2;
	real atp(time) micromolar;
	real weight2(time) dimensionless;
	real topa2 dimensionless;
	real bottom2(time) dimensionless;
	real topa3(time) dimensionless;
	real weight3(time) dimensionless;
	real bottom3(time) dimensionless;
	real famp dimensionless;
	famp=0.02;
	real fatp dimensionless;
	fatp=20;
	real ffbp dimensionless;
	ffbp=0.2;
	real fbt dimensionless;
	fbt=20;
	real fmt dimensionless;
	fmt=20;
	real weight4(time) dimensionless;
	real topa4(time) dimensionless;
	real bottom4(time) dimensionless;
	real weight5(time) dimensionless;
	real topa5(time) dimensionless;
	real bottom5(time) dimensionless;
	real weight6(time) dimensionless;
	real topa6(time) dimensionless;
	real bottom6(time) dimensionless;
	real weight7(time) dimensionless;
	real topa7(time) dimensionless;
	real bottom7(time) dimensionless;
	real weight8(time) dimensionless;
	real topa8(time) dimensionless;
	real bottom8(time) dimensionless;
	real weight9(time) dimensionless;
	real topa9(time) dimensionless;
	real bottom9(time) dimensionless;
	real weight10(time) dimensionless;
	real topa10(time) dimensionless;
	real bottom10(time) dimensionless;
	real weight11(time) dimensionless;
	real topa11(time) dimensionless;
	real bottom11(time) dimensionless;
	real weight12(time) dimensionless;
	real topa12(time) dimensionless;
	real bottom12(time) dimensionless;
	real weight13(time) dimensionless;
	real topa13(time) dimensionless;
	real bottom13(time) dimensionless;
	real weight14(time) dimensionless;
	real topa14(time) dimensionless;
	real bottom14(time) dimensionless;
	real weight15(time) dimensionless;
	real topa15(time) dimensionless;
	real bottom15(time) dimensionless;
	real weight16(time) dimensionless;
	real topa16(time) dimensionless;
	real bottom16(time) dimensionless;
	real topb(time) dimensionless;
	real amp(time) micromolar;
	real mgadp(time) micromolar;
	real adp3m(time) micromolar;
	real atp4m(time) micromolar;
	real topo(time) dimensionless;
	real bottomo(time) dimensionless;
	real katpo(time) dimensionless;
	real gkatpbar picoS;
	gkatpbar=25000;
	real adp(time) micromolar;
	when(time=time.min) adp=780;
	real kdd dimensionless;
	kdd=17;
	real ktd dimensionless;
	ktd=26;
	real ktt dimensionless;
	ktt=1;
	real fback(time) dimensionless;
	real taua dimensionless;
	taua=300000;
	real r1 micromolar;
	r1=0.35;
	real r dimensionless;
	r=1;
	real y(time) dimensionless;
	real vg dimensionless;
	vg=2.2;
	real kg flux;
	kg=10;
	real rad(time) dimensionless;
	real ratio(time) dimensionless;

	// <component name="environment">

	// <component name="membrane">
	V:time=((-1)*(IK+ICa+IKCa+IKATP)/Cm);
	IK=(gK*n*(V-VK));
	IKCa=(gKCa/(1+(kd/c)^2)*(V-VK));
	ICa=(gCa*minf*(V-VCa));
	minf=(1/(1+exp((-1)*(20+V/(1 millivolt))/12)));
	ninf=(1/(1+exp((-1)*(16+V/(1 millivolt))/5)));
	n:time=((ninf-n)/taun);

	// <component name="calcium_handling">
	c:time=(fcyt*(Jmem+Jer));
	cer:time=((-1)*fer*sigmaV*Jer);
	was_Jserca=(Kserca*c);
	Jleak=(pleak*(cer-c));
	Jer=(epser*(Jleak-was_Jserca)/lambdaer);
	Jmem=((-1)*(alpha*ICa+kpmca*c));

	// <component name="glycolysis">
	f6p=(.3*g6p);
	rgpdh=(.2*sqrt(abs(fbp*(1 micromolar)/(1 second)^2)));
	fbp:time=(lambda*(pfk/(1 second)-.5*rgpdh));
	g6p:time=(lambda*(Rgk*(1 micromolar)-pfk/(1 second)));

	// <component name="pfk">
	pfk=((1 micromolar)*((pfkbas*cat*topa16+cat*topb)/bottom16));
	weight2=(atp^2/(k4*(1 micromolar)));
	topa2=topa1;
	bottom2=(bottom1+weight2);
	weight3=(f6p^2/(k3*(1 micromolar)));
	topa3=(topa2+weight3);
	bottom3=(bottom2+weight3);
	weight4=((f6p*atp)^2/(fatp*k3*k4*(1 micromolar)^2));
	topa4=(topa3+weight4);
	bottom4=(bottom3+weight4);
	weight5=(fbp/k2);
	topa5=topa4;
	bottom5=(bottom4+weight5);
	weight6=(fbp*atp^2/(k2*k4*fbt*(1 micromolar)));
	topa6=topa5;
	bottom6=(bottom5+weight6);
	weight7=(fbp*f6p^2/(k2*k3*ffbp*(1 micromolar)));
	topa7=(topa6+weight7);
	bottom7=(bottom6+weight7);
	weight8=(fbp*f6p^2*atp^2/(k2*k3*k4*ffbp*fbt*fatp*(1 micromolar)^2));
	topa8=(topa7+weight8);
	bottom8=(bottom7+weight8);
	weight9=(amp/k1);
	topa9=topa8;
	bottom9=(bottom8+weight9);
	weight10=(amp*atp^2/(k1*k4*fmt*(1 micromolar)));
	topa10=topa9;
	bottom10=(bottom9+weight10);
	weight11=(amp*f6p^2/(k1*k3*famp*(1 micromolar)));
	topa11=(topa10+weight11);
	bottom11=(bottom10+weight11);
	weight12=(amp*f6p^2*atp^2/(k1*k3*k4*famp*fmt*fatp*(1 micromolar)^2));
	topa12=(topa11+weight12);
	bottom12=(bottom11+weight12);
	weight13=(amp*fbp/(k1*k2));
	topa13=topa12;
	bottom13=(bottom12+weight13);
	weight14=(amp*fbp*atp^2/(k1*k2*k4*fbt*fmt*(1 micromolar)));
	topa14=topa13;
	bottom14=(bottom13+weight14);
	weight15=(amp*fbp*f6p^2/(k1*k2*k3*ffbp*famp*(1 micromolar)));
	topa15=topa14;
	bottom15=(bottom14+weight15);
	topb=weight15;
	weight16=(amp*fbp*f6p^2*atp^2/(k1*k2*k3*k4*ffbp*famp*fbt*fmt*fatp*(1 micromolar)^2));
	topa16=(topa15+weight16);
	bottom16=(bottom15+weight16);

	// <component name="KATP">
	mgadp=(.165*adp);
	adp3m=(.135*adp);
	atp4m=(.05*atp);
	topo=(.08*(1+2*mgadp/(kdd*(1 micromolar)))+.89*(mgadp/(kdd*(1 micromolar)))^2);
	bottomo=((1+mgadp/(kdd*(1 micromolar)))^2*(1+adp3m/(ktd*(1 micromolar))+atp4m/(ktt*(1 micromolar))));
	katpo=(topo/bottomo);
	IKATP=(gkatpbar*katpo*(V-VK));

	// <component name="nucleotides">
	adp:time=((atp-adp*exp(fback*(1-c/r1)))/(taua*(1 millisecond)));
	fback=(r+y);
	y=(vg*(rgpdh/(kg+rgpdh)));
	atp=(.5*(atot-adp+rad*(1 micromolar)));
	amp=(adp*adp/atp);
	ratio=(atp/adp);
	rad=(sqrt(abs((adp-atot)^2-4*adp^2))/(1 micromolar));
}