import nsrunit;
unit conversion on;
unit s=1 second^1;
unit per_s=1 second^(-1);
unit uM=1E-3 meter^(-3)*mole^1;
unit uM_per_s=1E-3 meter^(-3)*second^(-1)*mole^1;

math main {
	realDomain time s;
	time.min=0;
	extern time.max;
	extern time.delta;
	real v_0 uM_per_s;
	v_0=1;
	real v_1 uM_per_s;
	v_1=7.3;
	real v_2(time) uM_per_s;
	real v_3(time) uM_per_s;
	real beta dimensionless;
	beta=0.301;
	real k_f per_s;
	k_f=1;
	real k per_s;
	k=10;
	real Y(time) uM;
	when(time=time.min) Y=1.75;
	real Z(time) uM;
	when(time=time.min) Z=0.5;
	real V_M2 uM_per_s;
	V_M2=65;
	real v_2.K_2 uM;
	v_2.K_2=1;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real n dimensionless;
//	n=2;
	real V_M3 uM_per_s;
	V_M3=500;
	real K_A uM;
	K_A=0.9;
	real K_R uM;
	K_R=2;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real m dimensionless;
//	m=2;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real p dimensionless;
//	p=4;
	real W_star(time) dimensionless;
	when(time=time.min) W_star=0;
	real W_T uM;
	W_T=1;
	real v_P uM_per_s;
	v_P=5;
	real v_K(time) uM_per_s;
	real K_1 dimensionless;
	K_1=0.1;
	real W_star.K_2 dimensionless;
	W_star.K_2=0.1;
	real V_MK uM_per_s;
	V_MK=40;
	real K_a uM;
	K_a=2.5;

	// <component name="environment">

	// <component name="Ca">
	Z:time=(v_0+v_1*beta-v_2+v_3+k_f*Y-k*Z);
	Y:time=(v_2-v_3-k_f*Y);

	// <component name="v_2">
	v_2=(V_M2*(Z^2/(v_2.K_2^2+Z^2)));

	// <component name="v_3">
	v_3=(V_M3*(Y^2/(K_R^2+Y^2))*(Z^4/(K_A^4+Z^4)));

	// <component name="W_star">
	W_star:time=(v_P/W_T*(v_K/v_P*((1-W_star)/(K_1+1+(-1)*W_star))-W_star/(W_star.K_2+W_star)));

	// <component name="v_K">
	v_K=(V_MK*(Z/(K_a+Z)));
}