/*
 * A phantom bursting mechanism for episodic bursting
 * 
 * Model Status
 * 
 * This CellML model runs in both OpenCell and COR to replicate
 * the published results (figure 4 in the original paper where
 * g_K1 = 20pS).
 * 
 * Model Structure
 * 
 * ABSTRACT: We describe a novel dynamic mechanism for episodic
 * or compound bursting oscillations, in which bursts of electrical
 * impulses are clustered together into episodes, separated by
 * long silent phases. We demonstrate the mechanism for episodic
 * bursting using a minimal mathematical model for "phantom bursting."
 * Depending on the location in parameter space, this model can
 * produce fast, medium, or slow bursting, or in the present case,
 * fast, slow, and episodic bursting. The episodic bursting is
 * modestly robust to noise and to parameter variation, and the
 * effect that noise has on the episodic bursting pattern is quite
 * different from that of an alternate episodic burst mechanism
 * in which the slow envelope is produced by metabolic oscillations.
 * This mechanism could account for episodic bursting produced
 * in endocrine cells or neurons, such as pancreatic islets or
 * gonadotropin releasing neurons of the hypothalamus.s.
 * 
 * The original paper reference is cited below.
 * 
 * A phantom bursting mechanism for episodic bursting, R. Bertram,
 * J. Rhoads and W.P. Cimbora, 2008, Bulletin of Mathematical Biology,
 * 70, 1979-1993. PubMed ID: 18648884
 * 
 * cell schematic for the model
 * 
 * [[Image file: bertram_2008.png]]
 * 
 * Schematic diagram of the cell model.
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 2 equation(s)
unit conversion off;
// unit millisecond predefined
// unit millivolt predefined
// unit femtofarad predefined
// unit picosiemens predefined
// unit femtoampere predefined

math main {
	realDomain time millisecond;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) millivolt;
	when(time=time.min) V=-40.0;
	real Cm femtofarad;
	Cm=4525.0;
	real Ica(time) femtoampere;
	real Is1(time) femtoampere;
	real Is2(time) femtoampere;
	real Il(time) femtoampere;
	real Ik(time) femtoampere;
	real gCa picosiemens;
	gCa=280.0;
	real VCa millivolt;
	VCa=100.0;
	real m_infinity(time) dimensionless;
	real vm millivolt;
	vm=-22.0;
	real sm millivolt;
	sm=7.5;
	real gs1 picosiemens;
	gs1=20.0;
	real VK millivolt;
	VK=-80.0;
	real s1(time) dimensionless;
	when(time=time.min) s1=0.9;
	real s1_infinity(time) dimensionless;
	real autos1 dimensionless;
	autos1=1;
	real s1knot dimensionless;
	s1knot=1;
	real tau_s1 millisecond;
	tau_s1=1000.0;
	real vs1 millivolt;
	vs1=-50.0;
	real ss1 millivolt;
	ss1=5;
	real gK picosiemens;
	gK=1300.0;
	real n(time) dimensionless;
	when(time=time.min) n=0.0;
	real n_infinity(time) dimensionless;
	real tau_n_bar millisecond;
	tau_n_bar=8.25;
	real tau_n(time) millisecond;
	real vn millivolt;
	vn=-9.0;
	real sn millivolt;
	sn=10.0;
	real gl picosiemens;
	gl=25.0;
	real Vl millivolt;
	Vl=-40.0;
	real gs2 picosiemens;
	gs2=16;
	real s2(time) dimensionless;
	when(time=time.min) s2=0.5;
	real s2_infinity(time) dimensionless;
	real autos2 dimensionless;
	autos2=1;
	real s2knot dimensionless;
	s2knot=0.49;
	real tau_s2 millisecond;
	tau_s2=30000.0;
	real vs2 millivolt;
	vs2=-40.0;
	real ss2 millivolt;
	ss2=15;

	// <component name="environment">

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

	// <component name="Ica">
	Ica=(gCa*m_infinity*(V-VCa));

	// <component name="m">
	m_infinity=(1/(1+exp((vm-V)/sm)));

	// <component name="Is1">
	Is1=(gs1*s1*(V-VK));

	// <component name="s1">
	s1:time=(autos1*((s1_infinity-s1)/tau_s1)+(1-autos1)*(s1knot-s1));
	s1_infinity=(1/(1+exp((vs1-V)/ss1)));

	// <component name="Ik">
	Ik=(gK*n*(V-VK));

	// <component name="n">
	n:time=((n_infinity-n)/tau_n);
	n_infinity=(1/(1+exp((vn-V)/sn)));
	tau_n=(tau_n_bar/(1+exp((V-vn)/sn)));

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

	// <component name="Is2">
	Is2=(gs2*s2*(V-VK));

	// <component name="s2">
	s2:time=(autos2*((s2_infinity-s2)/tau_s2)+(1-autos2)*(s2knot-s2));
	s2_infinity=(1/(1+exp((vs2-V)/ss2)));
}