/*
 * The Mechanical Properties of Cardiac Muscle, 1998
 * 
 * Model Status
 * 
 * This model is not currently functional (overconstrained).
 * 
 * Model Structure
 * 
 * Finite element models of the electrical and mechanical behaviour
 * of the whole heart are well established. At the cellular level,
 * there are several models published which describe the changes
 * in ion concentrations and membrane ionic currents underlying
 * the cardiac cell action potential. However, the Hunter-McCulloch-ter
 * Keurs 1998 model is the first to capture the mechanical properties
 * of actively contracting cardiac muscle. The model is based on
 * an extensive review of experimental data from a variety of preparations
 * and species. These experiments are interpreted with a four state
 * variable model which includes i) the passive elasticity of myocardial
 * tissue, ii) the rapid binding of Ca2+ to troponin C, iii) the
 * kinetics of tropomyosin movement and iv) the kinetics of crossbridge
 * tension development.
 * 
 * The complete original paper reference is cited below:
 * 
 * Modelling the mechanical properties of cardiac muscle, P.J.
 * Hunter, A.D. McCulloch and H.E.D.J. ter Keurs, 1998, Progress
 * in Biophysics and Molecular Biology, 69, 289-331. (A PDF version
 * of the article is available to Science Direct subscribers.)
 * PubMed ID: 9785944
 * 
 * The raw CellML description of the Hunter-McCulloch-ter Keurs
 * model can be downloaded in various formats as described in .
 * For an example of a more complete documentation for an electrophysiological
 * model, see The Hodgkin-Huxley Squid Axon Model, 1952.
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 4 equation(s)
unit conversion off;
// unit micromolar predefined
unit first_order_rate_constant=1 second^(-1);
unit second_order_rate_constant=1E3 meter^3*second^(-1)*mole^(-1);
// unit kilopascal predefined
unit per_kilopascal=.001 kilogram^(-1)*meter^1*second^2;

math main {
	//Warning:  the following variables were set 'extern' or given
	//  an initial value of '0' because the model would otherwise be
	//  underdetermined:  Ca_b, z, tau
	//Warning:  the following variables had initial values which were
	//  suppressed because the model would otherwise be overdetermined:
	//   n
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Ca_i(time) micromolar;
	when(time=time.min) Ca_i=10.0;
	real Ca_max micromolar;
	Ca_max=1.0;
	real tau_Ca second;
	tau_Ca=0.06;
	real Ca_o micromolar;
	Ca_o=0.01;
	real Ca_b(time) micromolar;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) Ca_b=0;
	real lambda dimensionless;
	lambda=2.6;
	real Ca_b_max micromolar;
	Ca_b_max=2.26;
	real rho_0 second_order_rate_constant;
	rho_0=100.0;
	real rho_1 first_order_rate_constant;
	rho_1=163.0;
	real To(time) kilopascal;
	real T(time) kilopascal;
	real z(time) dimensionless;
	//Warning:  Assuming zero initial condition; nothing provided in original CellML model.
	when(time=time.min) z=0;
	real C_50 micromolar;
	C_50=1.0;
	real pC_50 micromolar;
	real pC_50_ref micromolar;
	pC_50_ref=6.2;
	real n dimensionless;
	//Warning: CellML initial value suppressed to prevent overdetermining model.  Original initial value: n=4.5;
	real n_ref dimensionless;
	n_ref=6.9;
	real alpha_0 first_order_rate_constant;
	alpha_0=2.0;
	real T_ref kilopascal;
	T_ref=100.0;
	real beta_0 per_kilopascal;
	beta_0=1.45;
	real beta_1 per_kilopascal;
	beta_1=1.95;
	real beta_2 per_kilopascal;
	beta_2=0.31;
	real a dimensionless;
	a=0.5;
	real Q(time) dimensionless;
	real A1 dimensionless;
	A1=50.0;
	real A2 dimensionless;
	A2=175.0;
	real A3 dimensionless;
	A3=175.0;
	real alpha_1 first_order_rate_constant;
	alpha_1=33.0;
	real alpha_2 first_order_rate_constant;
	alpha_2=2850.0;
	real alpha_3 first_order_rate_constant;
	alpha_3=2850.0;
	extern real tau second;
	real dlambda_dt(time) first_order_rate_constant;

	// <component name="environment">

	// <component name="calcium_transient">
	Ca_i:time=(Ca_o+(Ca_max-Ca_o)*(time/tau_Ca)*exp((1-time)/tau_Ca));

	// <component name="TnC_Ca_binding_kinetics">
	Ca_b:time=(rho_0*Ca_i*(Ca_b_max-Ca_b)-rho_1*Ca_b*(1-T/(lambda*To)));

	// <component name="thin_filament_kinetics">
	z:time=(alpha_0*((Ca_b/C_50)^n*(1-z)-z));
	To=(T_ref*(1+beta_0*(lambda-1))*z);
	n=(n_ref*(1+beta_1*(lambda-1)));
	pC_50=(pC_50_ref*(1+beta_2*(lambda-1)));

	// <component name="crossbridge_kinetics">
	dlambda_dt=(alpha_1/A1*((T/To-1)/(T/To+a)));
	T=(To*((1+a*Q)/(1-Q)));
	Q=(A1*(exp((-1)*alpha_1*(time-tau))*dlambda_dt*tau)+A2*(exp((-1)*alpha_2*(time-tau))*dlambda_dt*tau)+A3*(exp((-1)*alpha_3*(time-tau))*dlambda_dt*tau));
}