/*
 * Comparison of simulated and measured calcium sparks in intact
 * skeletal muscle fibers of the frog
 * 
 * Model Status
 * 
 * This CellML model runs in OpenCell and COR and the units are
 * consistent. There are no simple figures of concentration against
 * time in the publication to compare the CellML model simulation
 * output against.
 * 
 * Model Structure
 * 
 * ABSTRACT: Calcium sparks in frog intact skeletal muscle fibers
 * were modeled as stereotypical events that arise from a constant
 * efflux of Ca(2+) from a point source for a fixed period of time
 * (e.g., 2.5 pA of Ca(2+) current for 4.6 ms; 18 degrees C). The
 * model calculates the local changes in the concentrations of
 * free Ca(2+) and of Ca(2+) bound to the major intrinsic myoplasmic
 * Ca(2+) buffers (troponin, ATP, parvalbumin, and the SR Ca(2+)
 * pump) and to the Ca(2+) indicator (fluo-3). A distinctive feature
 * of the model is the inclusion of a binding reaction between
 * fluo-3 and myoplasmic proteins, a process that strongly affects
 * fluo-3's Ca(2+)-reaction kinetics, its apparent diffusion constant,
 * and hence the morphology of sparks. DeltaF/F (the change in
 * fluo-3's fluorescence divided by its resting fluorescence) was
 * estimated from the calculated changes in fluo-3 convolved with
 * the microscope point-spread function. To facilitate comparisons
 * with measured sparks, noise and other sources of variability
 * were included in a random repetitive fashion to generate a large
 * number of simulated sparks that could be analyzed in the same
 * way as the measured sparks. In the initial simulations, the
 * binding of Ca(2+) to the two regulatory sites on troponin was
 * assumed to follow identical and independent binding reactions.
 * These simulations failed to accurately predict the falling phase
 * of the measured sparks. A second set of simulations, which incorporated
 * the idea of positive cooperativity in the binding of Ca(2+)
 * to troponin, produced reasonable agreement with the measurements.
 * Under the assumption that the single channel Ca(2+) current
 * of a ryanodine receptor (RYR) is 0.5-2 pA, the results suggest
 * that 1-5 active RYRs generate an average Ca(2+) spark in a frog
 * intact muscle fiber.
 * 
 * The original paper reference is cited below:
 * 
 * Comparison of simulated and measured calcium sparks in intact
 * skeletal muscle fibers of the frog, S.M. Baylor, S. Hollingworth
 * and W.K. Chandler, 2002, Journal of General Physiology, 120,
 * 349-368. PubMed ID: 12198091
 * 
 * [[Image file: baylor_2002.png]]
 * 
 * Schematic diagrams of the Ca2+ binding reactions for various
 * buffers and indicators: A The reaction of Ca2+ with ATP in the
 * presence of free Mg2+, B Reaction of Ca2+ with protein (Pr)
 * and fluo-3 (Fluo), C Competitive reaction of Ca2+ and Mg2+ with
 * parvalbumin (Parv), D Binding reaction of Ca2+ binding and transport
 * by the sarcoplasmic reticulum Ca2+ pump (E), E One-step reaction
 * of Ca2+ with Troponin (Trop), and F Two-step reaction of Ca2+
 * with Troponin (Trop).
 */

import nsrunit;
unit conversion on;
// unit micromolar predefined
unit first_order_rate_constant=1 second^(-1);
unit second_order_rate_constant=1E3 meter^3*second^(-1)*mole^(-1);

math main {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Ca(time) micromolar;
	when(time=time.min) Ca=0.05;
	real k1 second_order_rate_constant;
	k1=2.033E14;
	real k1_ first_order_rate_constant;
	k1_=2642.0;
	real k2 second_order_rate_constant;
	k2=1.017E14;
	real k2_ first_order_rate_constant;
	k2_=13.21;
	real Trop(time) micromolar;
	when(time=time.min) Trop=360.0;
	real CaTrop(time) micromolar;
	when(time=time.min) CaTrop=0.01;
	real Ca2Trop(time) micromolar;
	when(time=time.min) Ca2Trop=0.01;

	// <component name="environment">

	// <component name="Ca">
	Ca:time=(k1_*CaTrop+k2_*Ca2Trop-(k1*Ca*Trop+k2*Ca*CaTrop));

	// <component name="Trop">
	Trop:time=(k1_*CaTrop-k1*Ca*Trop);

	// <component name="CaTrop">
	CaTrop:time=(k1*Ca*Trop+k2_*Ca2Trop-(k1_*CaTrop+k2*Ca*CaTrop));

	// <component name="Ca2Trop">
	Ca2Trop:time=(k2*Ca*CaTrop-k2_*Ca2Trop);

	// <component name="reaction_constants">
}