/*
 * Mechanical regulation of cardiac muscle by coupling calcium
 * kinetics with cross-bridge cycling: a dynamic model
 * 
 * Model Status
 * 
 * This model does not run in OpenCell or COR because all variables
 * are not initialized and some variables are undefined. This is
 * only the case for the component force equations but it seems
 * that all the variables in that component cannot be calculated
 * because some of the equations in the original paper look like
 * partial differential equations which are beyond the scope of
 * CellML.
 * 
 * Model Structure
 * 
 * Abstract: This study describes the regulation of mechanical
 * activity in the intact cardiac muscle, the effects of the free
 * calcium transients and the mechanical constraints, and emphasizes
 * the central role of the troponin complex in regulating muscle
 * activity. A "loose coupling" between calcium binding to troponin
 * and cross-bridge cycling is stipulated, allowing the existence
 * of cross bridges in the strong conformation without having bound
 * calcium on the neighboring troponin. The model includes two
 * feedback mechanisms: 1) a positive feedback, or cooperativity,
 * in which the cycling cross bridges affect the affinity of troponin
 * for calcium, and 2) a negative mechanical feedback, where the
 * filament-sliding velocity affects cross-bridge cycling. The
 * model simulates the reported experimental force-length and force-velocity
 * relationships at different levels of activation. The dependence
 * of the shortening velocity on calcium concentration, sarcomere
 * length, internal load, and rate of cross-bridge cycling is described
 * analytically in agreement with reported data. Furthermore, the
 * model provides an analytic solution for Hill's equation of the
 * force-velocity relationship and for the phenomena of unloaded
 * shortening velocity and force deficit. The model-calculated
 * changes in free calcium in various mechanical conditions are
 * in good agreement with the available experimental results.
 * 
 * The complete original paper reference is cited below:
 * 
 * Mechanical regulation of cardiac muscle by coupling calcium
 * kinetics with cross-bridge cycling: a dynamic model, A. Landesberg,
 * S. Sideman, 1994, Am. J. Physiol- Heart Circ. Physiol, 267,
 * H779-795. PubMed ID: 8067434
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 18 equation(s)
unit conversion off;
unit umeter=1E-6 meter^1;
unit uM=1E-3 meter^(-3)*mole^1;
unit per_second=1 second^(-1);
unit uM_per_second=1E-3 meter^(-3)*second^(-1)*mole^1;
unit uM_per_umeter=1E3 meter^(-4)*mole^1;
unit uM_per_umeter_per_second=1E3 meter^(-4)*second^(-1)*mole^1;
unit per_uM_per_second=1E3 meter^3*second^(-1)*mole^(-1);
unit mN=.001 kilogram^1*meter^1*second^(-2);
unit umeter_per_second=1E-6 meter^1*second^(-1);
unit mN_second_per_umeter=1E3 kilogram^1*second^(-1);

math main {
	//Warning:  the following variables were set 'extern' or given
	//  an initial value of '0' because the model would otherwise be
	//  underdetermined:  F_bar, eta, eta_PE, force_equations.T_s, force_equations.U_s,
	//  E, D, B, SL, Sp_0
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real L_a umeter;
	L_a=1.15;
	real L_prime_m umeter;
	L_prime_m=1.5;
	real L_b umeter;
	L_b=0.1;
	real L_z umeter;
	L_z=0.1;
	real L_m umeter;
	real Troponin uM;
	Troponin=60;
	real k_h per_uM_per_second;
	k_h=1e8;
	real k_minus_h per_second;
	k_minus_h=0.33;
	real k_1 per_uM_per_second;
	k_1=4e7;
	real k_m per_uM_per_second;
	real K_1(time) per_second;
	real K_m(time) per_second;
	real K_max dimensionless;
	K_max=2e6;
	real f per_second;
	f=40;
	real g_0 per_second;
	g_0=12;
	real g_prime_0 per_second;
	real g_1 dimensionless;
	g_1=15;
	real g_prime_1 dimensionless;
	real k_minus_1 per_second;
	real k_minus_m per_second;
	real Ca(time) uM;
	when(time=time.min) Ca=0;
	real R_n(time) uM;
	real Rbar_n(time) uM_per_umeter;
	real A_n(time) uM;
	real Abar_n(time) uM_per_umeter;
	real alpha dimensionless;
	alpha=1;
	real R_s(time) uM;
	when(time=time.min) R_s=0;
	real Rbar_s(time) uM_per_umeter;
	when(time=time.min) Rbar_s=0;
	real dRbar_s_dt(time) uM_per_umeter_per_second;
	real A_s(time) uM;
	when(time=time.min) A_s=0;
	real Abar_s(time) uM_per_umeter;
	when(time=time.min) Abar_s=0;
	real dAbar_s_dt(time) uM_per_umeter_per_second;
	real single_overlap.T_s(time) uM;
	when(time=time.min) single_overlap.T_s=0;
	real Tbar_s(time) uM_per_umeter;
	when(time=time.min) Tbar_s=0;
	real dTbar_s_dt(time) uM_per_umeter_per_second;
	real single_overlap.U_s(time) uM;
	when(time=time.min) single_overlap.U_s=0;
	real Ubar_s(time) uM_per_umeter;
	when(time=time.min) Ubar_s=0;
	real dUbar_s_dt(time) uM_per_umeter_per_second;
	real Rbar_d(time) uM_per_umeter;
	when(time=time.min) Rbar_d=0;
	real Abar_d(time) uM_per_umeter;
	when(time=time.min) Abar_d=0;
	real Tbar_d(time) uM_per_umeter;
	when(time=time.min) Tbar_d=0;
	real Ubar_d(time) uM_per_umeter;
	when(time=time.min) Ubar_d=0;
	real V per_second;
	real R_d(time) uM;
	real A_d(time) uM;
	real T_d(time) uM;
	real U_d(time) uM;
	real Ca_0 uM;
	Ca_0=2000;
	real I_o(time) uM_per_second;
	real I_s(time) uM_per_second;
	real I_i(time) uM_per_second;
	real I_u(time) uM_per_second;
	real Q_o per_second;
	Q_o=30;
	real Q_s per_second;
	Q_s=0.3;
	real tau_SR second;
	tau_SR=0.002;
	real tau_SF second;
	tau_SF=0.02;
	real Q_i per_second;
	Q_i=0.6;
	real tau_iR second;
	tau_iR=0.002;
	real tau_iF second;
	tau_iF=0.1;
	real Q_u uM_per_second;
	Q_u=1000;
	real K_mu uM;
	K_mu=1;
	real I_l uM_per_second;
	I_l=0;
	real Ca_r uM;
	Ca_r=0;
	real BCa_L(time) uM;
	real BCa_h(time) uM;
	when(time=time.min) BCa_h=0;
	real F mN;
	real F_CB mN;
	real F_CE mN;
	real F_PE mN;
	extern real F_bar mN;
	extern real eta mN_second_per_umeter;
	extern real eta_PE mN_second_per_umeter;
	extern real force_equations.T_s dimensionless;
	extern real force_equations.U_s dimensionless;
	extern real E mN;
	extern real D dimensionless;
	extern real B mN;
	extern real SL umeter;
	extern real Sp_0 umeter;
	real d_alpha_dt per_second;
	d_alpha_dt=0;

	// <component name="environment">

	// <component name="parameters">
	L_m=(L_prime_m-L_b);
	K_1=(k_1*Ca);
	K_m=(k_m*Ca);
	k_m=k_1;
	k_minus_1=(k_1/K_max);
	k_minus_m=k_minus_1;
	g_prime_0=g_0;
	g_prime_1=g_1;

	// <component name="non_overlap">
	Rbar_n=(A_n*k_minus_1/(1+K_1));
	R_n=(Rbar_n*alpha*L_m);
	Abar_n=(R_n*K_1/(1+k_minus_1));
	A_n=(Abar_n*alpha*L_m);

	// <component name="single_overlap">
	dRbar_s_dt=((-1)*K_1*R_s+k_minus_1*A_s+(g_prime_0+g_prime_1*V)*single_overlap.U_s);
	Rbar_s:time=dRbar_s_dt;
	R_s:time=(if (V>(0 per_second)) dRbar_s_dt*alpha*L_m-Rbar_s*L_m*V else if (V<(0 per_second)) dRbar_s_dt*alpha*L_m+(2*Rbar_d-Rbar_s)*L_m*abs(V) else dRbar_s_dt*alpha*L_m);
	dAbar_s_dt=(K_1*R_s+((-1)*f-k_minus_1)*A_s+(g_0+g_1*V)*single_overlap.T_s);
	Abar_s:time=dAbar_s_dt;
	A_s:time=(if (V>(0 umeter_per_second)) dAbar_s_dt*alpha*L_m-Abar_s*L_m*V else if (V<(0 umeter_per_second)) dAbar_s_dt*alpha*L_m+(2*Abar_d-Abar_s)*L_m*abs(V) else dAbar_s_dt*alpha*L_m);
	dTbar_s_dt=(f*A_s+((-1)*g_0-g_1*V-k_minus_m)*single_overlap.T_s+K_m*single_overlap.U_s);
	Tbar_s:time=dTbar_s_dt;
	single_overlap.T_s:time=(if (V>(0 umeter_per_second)) dTbar_s_dt*alpha*L_m-Tbar_s*L_m*V else if (V<(0 umeter_per_second)) dTbar_s_dt*alpha*L_m+(2*Tbar_d-Tbar_s)*L_m*abs(V) else dTbar_s_dt*alpha*L_m);
	dUbar_s_dt=(k_minus_m*single_overlap.T_s+((-1)*K_m-g_prime_0-g_prime_1*V)*single_overlap.U_s);
	Ubar_s:time=dUbar_s_dt;
	single_overlap.U_s:time=(if (V>(0 umeter_per_second)) dUbar_s_dt*alpha*L_m-Ubar_s*L_m*V else if (V<(0 umeter_per_second)) dUbar_s_dt*alpha*L_m+(2*Ubar_d-Ubar_s)*L_m*abs(V) else dUbar_s_dt*alpha*L_m);

	// <component name="double_overlap">
	Rbar_d:time=((-1)*K_1*R_d+k_minus_1*A_d+(g_prime_0+g_prime_1*V)*U_d);
	R_d=(Rbar_d*alpha*L_m);
	Abar_d:time=(K_1*R_d+((-1)*f-k_minus_1)*A_d+(g_0+g_1*V)*T_d);
	A_d=(Abar_d*alpha*L_m);
	Tbar_d:time=(f*A_d+((-1)*g_0-g_1*V-k_minus_m)*T_d+K_m*U_d);
	T_d=(Tbar_d*alpha*L_m);
	Ubar_d:time=(k_minus_m*T_d+((-1)*K_m-g_prime_0-g_prime_1*V)*U_d);
	U_d=(Ubar_d*alpha*L_m);

	// <component name="calcium_flux">
	I_o=(Q_o*Ca);
	I_s=(Q_s*((1-exp((-1)*time/tau_SR))*exp((-1)*time/tau_SF)+I_l)*Ca_0);
	I_i=(Q_i*((1-exp((-1)*time/tau_iR))*exp((-1)*time/tau_iF)+I_l)*Ca_r);
	I_u=(Q_u*Ca/(K_mu+Ca));
	Ca:time=(I_s+I_i-(I_o+I_u));

	// <component name="bound_calcium">
	BCa_L=(A_s+single_overlap.T_s+A_d+T_d+A_n);
	BCa_h:time=((2*Troponin-BCa_h)*Ca*k_h-BCa_h*k_minus_h);

	// <component name="force_equations">
	F_CB=(F_bar-eta*V);
	F_CE=((force_equations.T_s+force_equations.U_s)*F_CB);
	F_PE=(if (SL>=Sp_0) E*(exp(D*(SL/Sp_0-1))-1)+eta_PE*V else (-1)*B*(1-SL/Sp_0)+eta_PE*V);
	F=(F_CE+F_PE);

	// <component name="input">
	V=d_alpha_dt;
}