/*
 * Voltage oscillations in the barnacle giant muscle fiber
 * 
 * Model Status
 * 
 * This CellML model represents the complete model (equations 1
 * and 2) with parameters values taken from the legend of figure
 * 6. The model runs in both OpenCell and COR, and the units are
 * consistent, however it does not replicate the published results
 * (figure 6).
 * 
 * Model Structure
 * 
 * ABSTRACT: Barnacle muscle fibers subjected to constant current
 * stimulation produce a variety of types of oscillatory behavior
 * when the internal medium contains the Ca++ chelator EGTA. Oscillations
 * are abolished if Ca++ is removed from the external medium, or
 * if the K+ conductance is blocked. Available voltage-clamp data
 * indicate that the cell's active conductance systems are exceptionally
 * simple. Given the complexity of barnacle fiber voltage behavior,
 * this seems paradoxical. This paper presents an analysis of the
 * possible modes of behavior available to a system of two noninactivating
 * conductance mechanisms, and indicates a good correspondence
 * to the types of behavior exhibited by barnacle fiber. The differential
 * equations of a simple equivalent circuit for the fiber are dealt
 * with by means of some of the mathematical techniques of nonlinear
 * mechanics. General features of the system are (a) a propensity
 * to produce damped or sustained oscillations over a rather broad
 * parameter range, and (b) considerable latitude in the shape
 * of the oscillatory potentials. It is concluded that for cells
 * subject to changeable parameters (either from cell to cell or
 * with time during cellular activity), a system dominated by two
 * noninactivating conductances can exhibit varied oscillatory
 * and bistable behavior.
 * 
 * The original paper reference is cited below:
 * 
 * Voltage oscillations in the barnacle giant muscle fiber, Morris
 * C and Lecar H, 1981, Biophysical Journal, 35, 193-213. PubMed
 * ID: 7260316
 * 
 * [[Image file: morris_1981.png]]
 * 
 * Schematic diagrams of the cell model.
 */

import nsrunit;
unit conversion on;
// unit millisecond predefined
unit per_millisecond=1E3 second^(-1);
// unit millivolt predefined
unit microF_per_cm2=.01 kilogram^(-1)*meter^(-4)*second^4*ampere^2;
unit microA_per_cm2=.01 meter^(-2)*ampere^1;
unit milliS_per_cm2=10 kilogram^(-1)*meter^(-4)*second^3*ampere^2;

math main {
	realDomain time millisecond;
	time.min=0;
	extern time.max;
	extern time.delta;
	real V(time) millivolt;
	when(time=time.min) V=-50;
	real C microF_per_cm2;
	C=20.0;
	real i_app microA_per_cm2;
	i_app=540.0;
	real i_L(time) microA_per_cm2;
	real i_Ca(time) microA_per_cm2;
	real i_K(time) microA_per_cm2;
	real g_L milliS_per_cm2;
	g_L=2.0;
	real E_L millivolt;
	E_L=-50.00;
	real E_Ca millivolt;
	E_Ca=100.0;
	real g_Ca milliS_per_cm2;
	g_Ca=4.0;
	real m(time) dimensionless;
	when(time=time.min) m=0.1;
	real m_infinity(time) dimensionless;
	real lambda_m(time) per_millisecond;
	real lambda_m_bar per_millisecond;
	lambda_m_bar=1.0;
	real V1 millivolt;
	V1=0.0;
	real V2 millivolt;
	V2=15.0;
	real E_K millivolt;
	E_K=-70.0;
	real g_K milliS_per_cm2;
	g_K=8.0;
	real n(time) dimensionless;
	when(time=time.min) n=0.1;
	real n_infinity(time) dimensionless;
	real lambda_n(time) per_millisecond;
	real lambda_n_bar per_millisecond;
	lambda_n_bar=0.1;
	real V3 millivolt;
	V3=10.0;
	real V4 millivolt;
	V4=10.0;

	// <component name="environment">

	// <component name="V">
	V:time=((i_app-(i_L+i_Ca+i_K))/C);

	// <component name="i_L">
	i_L=(g_L*(V-E_L));

	// <component name="i_Ca">
	i_Ca=(g_Ca*m*(V-E_Ca));

	// <component name="i_Ca_m_gate">
	m_infinity=(.5*(1+tanh((V-V1)/V2)));
	lambda_m=(lambda_m_bar*cosh((V-V1)/(2*V2)));
	m:time=(lambda_m*(m_infinity-m));

	// <component name="i_K">
	i_K=(g_K*n*(V-E_K));

	// <component name="i_K_n_gate">
	n_infinity=(.5*(1+tanh((V-V3)/V4)));
	lambda_n=(lambda_n_bar*cosh((V-V3)/(2*V4)));
	n:time=(lambda_n*(n_infinity-n));
}