/*
 * The phantom burster model for pancreatic beta-cells
 * 
 * Model Status
 * 
 * This model has been rebuilt according to the author's original
 * XPPAUT code, which can be found here. This version of the CellML
 * model represents the fast bursting model where gs1=20. The model
 * replicates figure 2 in the published paper. The model runs in
 * both PCEnv and COR and the units are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: Pancreatic beta-cells exhibit bursting oscillations
 * with a wide range of periods. Whereas periods in isolated cells
 * are generally either a few seconds or a few minutes, in intact
 * islets of Langerhans they are intermediate (10-60 s). We develop
 * a mathematical model for beta-cell electrical activity capable
 * of generating this wide range of bursting oscillations. Unlike
 * previous models, bursting is driven by the interaction of two
 * slow processes, one with a relatively small time constant (1-5
 * s) and the other with a much larger time constant (1-2 min).
 * Bursting on the intermediate time scale is generated without
 * need for a slow process having an intermediate time constant,
 * hence phantom bursting. The model suggests that isolated cells
 * exhibiting a fast pattern may nonetheless possess slower processes
 * that can be brought out by injecting suitable exogenous currents.
 * Guided by this, we devise an experimental protocol using the
 * dynamic clamp technique that reliably elicits islet-like, medium
 * period oscillations from isolated cells. Finally, we show that
 * strong electrical coupling between a fast burster and a slow
 * burster can produce synchronized medium bursting, suggesting
 * that islets may be composed of cells that are intrinsically
 * either fast or slow, with few or none that are intrinsically
 * medium.
 * 
 * The original paper reference is cited below:
 * 
 * The phantom burster model for pancreatic beta-cells, Richard
 * Bertram, Joseph Previte, Arthur Sherman, Tracie A. Kinard and
 * Leslie S. Satin, 2000, Biophysical Journal, 79, 2880-2892. PubMed
 * ID: 11106596
 * 
 * cell schematic for the model
 * 
 * [[Image file: bertram_2000.png]]
 * 
 * Schematic diagram of the pancreatic beta-cell plasma membrane
 * showing the ionic currents captured by the phantom burster model.
 */

import nsrunit;
unit conversion on;
unit minute=60 second^1;
unit femtoA=1E-15 ampere^1;
unit femtoF=1E-15 kilogram^(-1)*meter^(-2)*second^4*ampere^2;
// unit millivolt predefined
unit picoS=1E-12 kilogram^(-1)*meter^(-2)*second^3*ampere^2;
// unit millisecond predefined

math main {
	realDomain time millisecond;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Cm femtoF;
	Cm=4524;
	real V(time) millivolt;
	when(time=time.min) V=-43;
	real ICa(time) femtoA;
	real IK(time) femtoA;
	real Il(time) femtoA;
	real Is1(time) femtoA;
	real Is2(time) femtoA;
	real Vm millivolt;
	Vm=-22;
	real VCa millivolt;
	VCa=100;
	real gCa picoS;
	gCa=280;
	real minf(time) dimensionless;
	real sm millivolt;
	sm=7.5;
	real VK millivolt;
	VK=-80;
	real gK picoS;
	gK=1300;
	real n(time) dimensionless;
	when(time=time.min) n=0.03;
	real lambda dimensionless;
	lambda=1.1;
	real tnbar dimensionless;
	tnbar=9.09;
	real Vn millivolt;
	Vn=-9;
	real sn millivolt;
	sn=10;
	real taun(time) dimensionless;
	real ninf(time) dimensionless;
	real gs1 picoS;
	gs1=20;
	real s1(time) dimensionless;
	when(time=time.min) s1=0.1;
	real s1inf(time) dimensionless;
	real Vs1 millivolt;
	Vs1=-40;
	real ss1 millivolt;
	ss1=0.5;
	real taus1 dimensionless;
	taus1=1000;
	real Vs2 millivolt;
	Vs2=-42;
	real s2(time) dimensionless;
	when(time=time.min) s2=0.434;
	real s2inf(time) dimensionless;
	real ss2 millivolt;
	ss2=0.4;
	real gs2 picoS;
	gs2=32;
	real taus2 dimensionless;
	taus2=120000;
	real gl picoS;
	gl=25;
	real Vl millivolt;
	Vl=-40;

	// <component name="environment">

	// <component name="membrane">
	V:time=((-1)*(ICa+IK+Il+Is1+Is2)/Cm);

	// <component name="Ca_current">
	minf=(1/(1+exp((Vm-V)/sm)));
	ICa=(gCa*minf*(V-VCa));

	// <component name="rapid_K_current">
	IK=(gK*n*(V-VK));
	n:time=(lambda*(ninf-n)/(taun*(1 millisecond)));
	ninf=(1/(1+exp((Vn-V)/sn)));
	taun=(tnbar/(1+exp((V-Vn)/sn)));

	// <component name="slow_K_current">
	Is1=(gs1*s1*(V-VK));
	s1inf=(1/(1+exp((Vs1-V)/ss1)));
	s1:time=((s1inf-s1)/(taus1*(1 millisecond)));

	// <component name="very_slow_K_current">
	Is2=(gs2*s2*(V-VK));
	s2inf=(1/(1+exp((Vs2-V)/ss2)));
	s2:time=((s2inf-s2)/(taus2*(1 millisecond)));

	// <component name="leak_current">
	Il=(gl*(V-Vl));
}