/*
 * A Bifurcation Analysis of Two Coupled Calcium Oscillators
 * 
 * Model Status
 * 
 * This CellML model runs in both OpenCell and COR to produce an
 * oscillating output similar to that from the original published
 * model. The units have been checked and they are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: In many cell types, asynchronous or synchronous oscillations
 * in the concentration of intracellular free calcium occur in
 * adjacent cells that are coupled by gap junctions. Such oscillations
 * are believed to underlie oscillatory intercellular calcium waves
 * in some cell types, and thus it is important to understand how
 * they occur and are modified by intercellular coupling. Using
 * a previous model of intracellular calcium oscillations in pancreatic
 * acinar cells, this article explores the effects of coupling
 * two cells with a simple linear diffusion term. Depending on
 * the concentration of a signal molecule, inositol (1,4,5)-trisphosphate,
 * coupling two identical cells by diffusion can give rise to synchronized
 * in-phase oscillations, as well as different-amplitude in-phase
 * oscillations and same-amplitude antiphase oscillations. Coupling
 * two nonidentical cells leads to more complex behaviors such
 * as cascades of period doubling and multiply periodic solutions.
 * This study is a first step towards understanding the role and
 * significance of the diffusion of calcium through gap junctions
 * in the coordination of oscillatory calcium waves in a variety
 * of cell types.
 * 
 * The original paper reference is cited below:
 * 
 * A bifurcation analysis of two coupled calcium oscillators, Michael
 * Bindschadler and James Sneyd, 2001, CHAOS, 11, 237-246. PubMed
 * ID: 12779457
 * 
 * cell diagram
 * 
 * [[Image file: bindschadler_2001.png]]
 * 
 * Schematic diagram of the IP3 receptor model. The receptor with
 * its three possible states: X, Y, and Z; representing open, shut
 * and inactive respectively, is embedded within a model of intracellular
 * calcium dynamics.
 */

import nsrunit;
unit conversion on;
unit per_second=1 second^(-1);
unit micro_molar=1E-3 meter^(-3)*mole^1;
unit micro_molar_per_second=1E-3 meter^(-3)*second^(-1)*mole^1;
unit second_order_rate=1E3 meter^3*second^(-1)*mole^(-1);

math main {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real phi3_c1(time) per_second;
	real h1(time) dimensionless;
	when(time=time.min) h1=0.8;
	real phi1_c1(time) second_order_rate;
	real phi2_c1(time) per_second;
	real p micro_molar;
	p=0.2778;
	real phi_1_c1(time) per_second;
	real phi3_c2(time) per_second;
	real h2(time) dimensionless;
	when(time=time.min) h2=0.1;
	real phi1_c2(time) second_order_rate;
	real phi2_c2(time) per_second;
	real phi_1_c2(time) per_second;
	real r2 second_order_rate;
	r2=100;
	real R1 micro_molar;
	R1=6;
	real k1 micro_molar_per_second;
	k1=44;
	real R3 micro_molar;
	R3=50;
	real k2 micro_molar_per_second;
	k2=26.5;
	real r4 per_second;
	r4=20;
	real k3 micro_molar_per_second;
	k3=1.6;
	real R5 micro_molar;
	R5=1.6;
	real c1(time) micro_molar;
	when(time=time.min) c1=0.3;
	real c2(time) micro_molar;
	when(time=time.min) c2=0.1;
	real Vp micro_molar_per_second;
	Vp=1.2;
	real Kp micro_molar;
	Kp=0.18;
	real j_pump_c1(time) micro_molar_per_second;
	real j_pump_c2(time) micro_molar_per_second;
	real kf micro_molar_per_second;
	kf=28;
	real j_receptor_c1(time) micro_molar_per_second;
	real j_receptor_c2(time) micro_molar_per_second;
	real j_diffusion(time) micro_molar_per_second;
	real D per_second;
	D=0.01;
	real j_leak micro_molar_per_second;
	j_leak=0.2;

	// <component name="environment">

	// <component name="h1">
	h1:time=(phi3_c1*(1-h1)-phi1_c1*phi2_c1*h1*p/(phi1_c1*p+phi_1_c1));

	// <component name="h2">
	h2:time=(phi3_c2*(1-h2)-phi1_c2*phi2_c2*h2*p/(phi1_c2*p+phi_1_c2));

	// <component name="phi">
	phi1_c1=(r2*c1/(R1+c1));
	phi_1_c1=(k1/(R3+c1));
	phi2_c1=((k2+r4*c1)/(R3+c1));
	phi3_c1=(k3/(R5+c1));
	phi1_c2=(r2*c2/(R1+c2));
	phi_1_c2=(k1/(R3+c2));
	phi2_c2=((k2+r4*c2)/(R3+c2));
	phi3_c2=(k3/(R5+c2));

	// <component name="j_pump">
	j_pump_c1=(Vp*c1^2/(Kp^2+c1^2));
	j_pump_c2=(Vp*c2^2/(Kp^2+c2^2));

	// <component name="j_receptor">
	j_receptor_c1=(kf*(p*h1*phi1_c1/(phi1_c1*p+phi_1_c1))^4);
	j_receptor_c2=(kf*(p*h2*phi1_c2/(phi1_c2*p+phi_1_c2))^4);

	// <component name="j_diffusion">
	j_diffusion=(D*(c2-c1));

	// <component name="c1">
	c1:time=(j_receptor_c1-j_pump_c1+j_leak+j_diffusion);

	// <component name="c2">
	c2:time=(j_receptor_c2-j_pump_c2+j_leak+j_diffusion);

	// <component name="model_parameters">
}