/*
 * Influence of viscosity on myocardium mechanical activity: a
 * mathematical model
 * 
 * Model Status
 * 
 * This model has been checked and it runs in both COR and PCEnv
 * to recreate the published results. Units are consistent throughout.
 * In this particular version of the model, time is expressed in
 * seconds and all calcium variables and parameters have been normalized
 * to the total concentration of TnC (parameter A_tot).
 * 
 * With default setting the model produce isometric contraction
 * (parameter isotonic=0). To switch from isomentric mode to isotonic
 * mode isotonic parameter should be set 1. Then set parameter
 * F_afterload to the value less then the maximum force calculated
 * at isometric mode. Make sure that the model is run in isometric
 * mode first and see the value of the maximum of isometric force.
 * If F_afterload is not less than this value, model performs isometric
 * mode instead of isotonic.
 * 
 * Model Structure
 * 
 * Abstract: We have previously proposed and validated a mathematical
 * model of myocardium contraction-relaxation cycle based on current
 * knowledge of regulatory role of Ca2+ and cross-bridge kinetics
 * in cardiac cell. That model did not include viscous elements.
 * Here we propose a modification of the model, in which two viscous
 * elements are added, one in parallel to the contractile element,
 * and one more in parallel to the series elastic element. The
 * modified model allowed us to simulate and explain some subtle
 * experimental data on relaxation velocity in isotonic twitches
 * and on a mismatch between the time course of sarcomere shortening/lengthening
 * and the time course of active force generation in isometric
 * twitches. Model results were compared with experimental data
 * obtained from 28 rat LV papillary muscles contracting and relaxing
 * against various loads. Additional model analysis suggested contribution
 * of viscosity to main inotropic and lusitropic characteristics
 * of myocardium performance.
 * 
 * model diagram Schematic diagram of the Katsnelson et al model
 * - an updated rheological scheme of the heart muscle including
 * contractile element CE, two passive elastic elements PE (parallel
 * one) and SE (series one) and two viscous elements VS1 and VS2.
 * The effect of calcium and the calcium binding ligand (B) in
 * facilitating actin-myosin binding is also highlighted.
 * 
 * The original paper reference is cited below:
 * 
 * Influence of viscosity on myocardium mechanical activity: a
 * mathematical model, Leonid B. Katsnelson, Larissa V. Nikitina,
 * Denis Chemla, Olga Solovyova, Catherine Coirault, Yves Lecarpentier
 * and Vladimir S. Markhasin, 2004, Journal of Theoretical Biology,
 * 230, 385-485. PubMed ID: 15302547
 */

import nsrunit;
unit conversion on;
unit per_second=1 second^(-1);
unit per_second2=1 second^(-2);
// unit micrometre predefined
unit micrometre_per_second=1E-6 meter^1*second^(-1);
unit micrometre_per_second2=1E-6 meter^1*second^(-2);
unit per_micrometre=1E6 meter^(-1);
unit second_per_micrometre=1E6 meter^(-1)*second^1;
// unit millinewton predefined
unit millinewton_per_second=.001 kilogram^1*meter^1*second^(-3);
unit millinewton_second_per_micrometre=1E3 kilogram^1*second^(-1);

math main {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real isotonic dimensionless;
	isotonic=0;
	real alpha_1 per_micrometre;
	alpha_1=19;
	real beta_1 millinewton;
	beta_1=0.29;
	real alpha_2 per_micrometre;
	alpha_2=14.6;
	real beta_2 millinewton;
	beta_2=0.000924;
	real alpha_3 per_micrometre;
	alpha_3=48;
	real beta_3 millinewton;
	beta_3=0.01;
	real lambda millinewton;
	lambda=96;
	real A_half dimensionless;
	A_half=0.6;
	real mu dimensionless;
	mu=3;
	real chi dimensionless;
	chi=0.705;
	real chi_0 dimensionless;
	chi_0=3;
	real m_0 dimensionless;
	m_0=0.9;
	real v_max micrometre_per_second;
	v_max=5.6;
	real a dimensionless;
	a=0.25;
	real d_h dimensionless;
	d_h=0.5;
	real alpha_P dimensionless;
	alpha_P=4;
	real l(time) micrometre;
	real F_muscle(time) millinewton;
	real flag(time) dimensionless;
	real F_afterload millinewton;
	F_afterload=2;
	real isotonic_mode(time) dimensionless;
	real l_0 micrometre;
	l_0=0.527;
	real S_0 micrometre;
	S_0=1.14;
	real q_v(time) per_second;
	real q_1 per_second;
	q_1=17.3;
	real q_2 per_second;
	q_2=259;
	real q_3 per_second;
	q_3=17.3;
	real q_4 per_second;
	q_4=15;
	real v_star micrometre_per_second;
	v_star=5.3035675;
	real v_1 micrometre_per_second;
	real alpha_G dimensionless;
	alpha_G=1;
	real a_on per_second;
	a_on=2300;
	real a_off per_second;
	a_off=290;
	real k_A dimensionless;
	k_A=2.8;
	real v(time) micrometre_per_second;
	when(time=time.min) v=0;
	real alpha_Q dimensionless;
	alpha_Q=10;
	real beta_Q dimensionless;
	beta_Q=5000;
	real F_CE(time) millinewton;
	real F_XSE(time) millinewton;
	real F_SE(time) millinewton;
	real F_PE(time) millinewton;
	real N(time) dimensionless;
	when(time=time.min) N=0.0001;
	real k_P_vis(time) millinewton_second_per_micrometre;
	real k_S_vis(time) millinewton_second_per_micrometre;
	real w(time) micrometre_per_second;
	when(time=time.min) w=0;
	real l_1(time) micrometre;
	when(time=time.min) l_1=0.437;
	real l_2(time) micrometre;
	when(time=time.min) l_2=0.439;
	real l_3(time) micrometre;
	when(time=time.min) l_3=0.089;
	real p_v(time) dimensionless;
	real K_chi(time) per_second;
	real M_A(time) dimensionless;
	real n_1(time) dimensionless;
	real L_oz(time) dimensionless;
	real k_p_v(time) per_second;
	real k_m_v(time) per_second;
	real A(time) dimensionless;
	when(time=time.min) A=0.01;
	real G_star(time) dimensionless;
	real P_star(time) dimensionless;
	real v_0 micrometre_per_second;
	v_0=560;
	real q_star per_second;
	q_star=1000;
	real dl_1_dt(time) micrometre_per_second;
	real dl_2_dt(time) micrometre_per_second;
	real dl_3_dt(time) micrometre_per_second;
	real phi_chi_2(time) micrometre_per_second;
	real phi_chi(time) micrometre_per_second2;
	real p_prime_v(time) second_per_micrometre;
	real alpha_P_lengthening per_micrometre;
	alpha_P_lengthening=16;
	real beta_P_lengthening millinewton_second_per_micrometre;
	beta_P_lengthening=0.0015;
	real alpha_P_shortening per_micrometre;
	alpha_P_shortening=16;
	real beta_P_shortening millinewton_second_per_micrometre;
	beta_P_shortening=0.0015;
	real alp_p(time) per_micrometre;
	real alpha_S_lengthening per_micrometre;
	alpha_S_lengthening=39;
	real beta_S_lengthening millinewton_second_per_micrometre;
	beta_S_lengthening=0.008;
	real alpha_S_shortening per_micrometre;
	alpha_S_shortening=46;
	real beta_S_shortening millinewton_second_per_micrometre;
	beta_S_shortening=0.006;
	real alp_s(time) per_micrometre;
	real gamma dimensionless;
	real case_1 second_per_micrometre;
	real case_2(time) second_per_micrometre;
	real case_3 second_per_micrometre;
	real case_4(time) second_per_micrometre;
	real dA_dt(time) per_second;
	real N_A(time) dimensionless;
	real pi_N_A(time) dimensionless;
	real B(time) dimensionless;
	when(time=time.min) B=0;
	real dB_dt(time) per_second;
	real Ca_C(time) dimensionless;
	when(time=time.min) Ca_C=0;
	real A_tot dimensionless;
	A_tot=1;
	real B_tot dimensionless;
	B_tot=0.4;
	real b_on per_second;
	b_on=2600;
	real b_off per_second;
	b_off=182;
	real a_c per_second2;
	a_c=5200;
	real r_Ca per_second;
	r_Ca=650;
	real q_Ca dimensionless;
	q_Ca=50;
	real t_d second;
	t_d=0.033;
	real Ca_m dimensionless;
	Ca_m=0.03;

	// <component name="environment">

	// <component name="parameters">

	// <component name="isotonic">
	flag=(if ((l>=l_0) and (time<(.15 second))) 0 else 1);
	isotonic_mode=(if (isotonic=0) 0 else if ((isotonic=1) and (F_muscle>=F_afterload)) 1 else if (((isotonic=1) and (l>=l_0)) and (flag=1)) 0 else 0);

	// <component name="parameters_izakov_et_al_1991">
	q_v=(if (v<=(0 micrometre_per_second)) q_1-q_2*v/v_max else if ((v<=v_star) and ((0 micrometre_per_second)<v)) (q_4-q_3)*v/v_star+q_3 else q_4/(1+beta_Q*(v-v_star)/v_max)^alpha_Q);
	v_1=(v_max/10);

	// <component name="force">
	F_CE=(lambda*p_v*N);
	F_muscle=F_XSE;
	F_SE=(beta_1*(exp(alpha_1*(l_2-l_1))-1));
	F_PE=(beta_2*(exp(alpha_2*l_2)-1));
	F_XSE=(beta_3*(exp(alpha_3*l_3)-1));

	// <component name="crossbridge_kinetics">
	M_A=(A^mu/(A^mu+A_half^mu));
	n_1=((.6 per_micrometre)*l_1+.5);
	L_oz=((l_1+S_0)/((.46 micrometre)+S_0));
	k_p_v=(chi*chi_0*q_v*m_0*G_star);
	k_m_v=(if (v<=v_star) chi_0*q_v*(1-chi*m_0*G_star) else chi_0*(q_4*(1-chi*m_0*G_star)+q_star*(v-v_star)/(v_0-v_star)));
	K_chi=(k_p_v*M_A*n_1*L_oz*(1-N)-k_m_v*N);
	N:time=K_chi;

	// <component name="length">
	l=(l_2+l_3);
	dl_1_dt=v;
	l_1:time=dl_1_dt;
	dl_2_dt=(if (k_S_vis=(0 millinewton_second_per_micrometre)) phi_chi_2 else w);
	l_2:time=dl_2_dt;
	dl_3_dt=(if (isotonic_mode=1) (0 micrometre_per_second) else if ((isotonic_mode=0) and (k_S_vis=(0 millinewton_second_per_micrometre))) (-1)*phi_chi_2 else (-1)*w);
	l_3:time=dl_3_dt;

	// <component name="CE_velocity">
	alp_p=(if (v<=(0 micrometre_per_second)) alpha_P_lengthening else alpha_P_shortening);
	k_P_vis=(if (v<=(0 micrometre_per_second)) beta_P_lengthening*exp(alpha_P_lengthening*l_1) else beta_P_shortening*exp(alpha_P_shortening*l_1));
	phi_chi=(if (isotonic_mode=1) (-1)*(lambda*K_chi*p_v+alp_p*k_P_vis*v^2+alpha_2*beta_2*exp(alpha_2*l_2)*w)/(lambda*N*p_prime_v+k_P_vis) else (-1)*(lambda*K_chi*p_v+alp_p*k_P_vis*v^2+(alpha_2*beta_2*exp(alpha_2*l_2)+alpha_3*beta_3*exp(alpha_3*l_3))*w)/(lambda*N*p_prime_v+k_P_vis));
	phi_chi_2=(if (isotonic_mode=1) alpha_1*beta_1*exp(alpha_1*(l_2-l_1))*v/(alpha_1*beta_1*exp(alpha_1*(l_2-l_1))+alpha_2*beta_2*exp(alpha_2*l_2)) else alpha_1*beta_1*exp(alpha_1*(l_2-l_1))*v/(alpha_1*beta_1*exp(alpha_1*(l_2-l_1))+alpha_2*beta_2*exp(alpha_2*l_2)+alpha_3*beta_3*exp(alpha_3*l_3)));
	v:time=(if (k_S_vis=(0 millinewton_second_per_micrometre)) (alpha_1*beta_1*exp(alpha_1*(l_2-l_1))*(phi_chi_2-v)-(lambda*K_chi*p_v+alp_p*k_P_vis*v^2))/(lambda*N*p_prime_v+k_P_vis) else phi_chi);

	// <component name="PE_velocity">
	alp_s=(if (w<=v) alpha_S_lengthening else alpha_S_shortening);
	k_S_vis=(if (w<=v) beta_S_lengthening*exp(alpha_S_lengthening*(l_2-l_1)) else beta_S_shortening*exp(alpha_S_shortening*(l_2-l_1)));
	w:time=(if (isotonic_mode=1) (k_S_vis*(phi_chi-alp_s*(w-v)^2)-alpha_1*beta_1*exp(alpha_1*(l_2-l_1))*(w-v)-alpha_2*beta_2*exp(alpha_2*l_2)*w)/k_S_vis else phi_chi-alp_s*(w-v)^2-(alpha_1*beta_1*exp(alpha_1*(l_2-l_1))*(w-v)+(alpha_2*beta_2*exp(alpha_2*l_2)+alpha_3*beta_3*exp(alpha_3*l_3))*w)/k_S_vis);

	// <component name="average_crossbridge_force">
	gamma=(a*d_h*(v_1/v_max)^2/(3*a*d_h-(a+1)*v_1/v_max));
	P_star=(if (v<=(0 micrometre_per_second)) a*(1+v/v_max)/(a-v/v_max) else 1+d_h-d_h^2*a/(a*d_h/gamma*(v/v_max)^2+(a+1)*v/v_max+a*d_h));
	G_star=(if ((((-1)*v_max)<=v) and (v<=(0 micrometre_per_second))) 1+.6*v/v_max else if (((0 micrometre_per_second)<v) and (v<=v_1)) P_star/((.4*a+1)*v/(a*v_max)+1) else P_star*exp((-1)*alpha_G*((v-v_1)/v_max)^alpha_P)/((.4*a+1)*v/(a*v_max)+1));
	case_1=(a*(.4+.4*a)/(v_max*((a+1)*.4)^2));
	case_2=(a*1*(1+.4*a+1.2*v/v_max+.6*(v/v_max)^2)/(v_max*((a-v/v_max)*(1+.6*v/v_max))^2));
	case_3=((.4*a+1)/(a*v_max));
	case_4=(1/v_max*exp((-1)*alpha_G*(v/v_max-v_1/v_max)^alpha_P)*((.4*a+1)/a+alpha_G*alpha_P*(1+(.4*a+1)*v/(a*v_max))*(v/v_max-v_1/v_max)^(alpha_P-1)));
	p_v=(P_star/G_star);
	p_prime_v=(if (v<=((-1)*v_max)) case_1 else if ((((-1)*v_max)<v) and (v<=(0 micrometre_per_second))) case_2 else if (((0 micrometre_per_second)<v) and (v<=v_1)) case_3 else case_4);

	// <component name="calcium_handling">
	N_A=(N/(L_oz*A));
	pi_N_A=(if (N_A>=1) 1 else .02^N_A);
	dA_dt=(a_on*(A_tot-A)*Ca_C-a_off*exp((-1)*k_A*A)*pi_N_A*A);
	A:time=dA_dt;
	dB_dt=(b_on*(B_tot-B)*Ca_C-b_off*B);
	B:time=dB_dt;
	Ca_C:time=(if (time<t_d) 4*a_c*Ca_m*time*(1-exp((-1)*a_c*time^2))*exp((-1)*a_c*time^2) else (-1)*dA_dt-dB_dt-r_Ca*exp((-1)*q_Ca*Ca_C)*Ca_C);
}