import nsrunit;
unit conversion on;
// unit nanomolar predefined
unit minute=60 second^1;
unit flux=1.6666667E-8 meter^(-3)*second^(-1)*mole^1;
unit first_order_rate_constant=.01666667 second^(-1);
unit second_order_rate_constant=1.6666667E4 meter^3*second^(-1)*mole^(-1);

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real RA(time) nanomolar;
	when(time=time.min) RA=0.1;
	real v_s1 flux;
	real k_d1 second_order_rate_constant;
	k_d1=1;
	real RA.C nanomolar;
	RA.C=0.1;
	real k_d5 first_order_rate_constant;
	k_d5=0;
	real M_C(time) nanomolar;
	when(time=time.min) M_C=0.1;
	real V_0 flux;
	V_0=0.365;
	real V_sC flux;
	V_sC=7.1;
	real F(time) nanomolar;
	when(time=time.min) F=0.0001;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real n dimensionless;
//	n=2;
	real K_A nanomolar;
	K_A=0.2;
	real k_d3 first_order_rate_constant;
	k_d3=1;
	real C.C(time) nanomolar;
	when(time=time.min) C.C=0.1;
	real k_s2 first_order_rate_constant;
	k_s2=1;
	real k_d2 first_order_rate_constant;
	k_d2=0.28;
	real k_s3 first_order_rate_constant;
	k_s3=1;
	real M_F nanomolar;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real m dimensionless;
//	m=2;
	real K_I nanomolar;
	K_I=0.2;
	real k_d4 first_order_rate_constant;
	k_d4=1;
	real k_s1 first_order_rate_constant;
	k_s1=1;
	real RALDH2_0 nanomolar;
	RALDH2_0=7.1;
	real x dimensionless;
	x=15;
	real L dimensionless;
	L=50;
	real M_0 nanomolar;
	M_0=5;
	real alpha_1(time) dimensionless;
	real K_r1 nanomolar;
	K_r1=1;
	real alpha_2(time) dimensionless;
	real K_r2 nanomolar;
	K_r2=1;
	real rho(time) dimensionless;

	// <component name="environment">

	// <component name="RA">
	RA:time=(v_s1-k_d1*RA.C*RA-k_d5*RA);

	// <component name="M_C">
	M_C:time=(V_0+V_sC*F^2/(K_A^2+F^2)-k_d3*M_C);

	// <component name="C">
	C.C:time=(k_s2*M_C-k_d2*C.C);

	// <component name="F">
	F:time=(k_s3*M_F*K_I^2/(K_I^2+RA^2)-k_d4*F);

	// <component name="v_s1">
	v_s1=(k_s1*RALDH2_0*(1-x/L));

	// <component name="M_F">
	M_F=(M_0*x/L);

	// <component name="alpha_1">

	// <component name="alpha_2">

	// <component name="rho">
	rho=(alpha_2/alpha_1);

	// <component name="parameters">

	// <component name="law1_alpha">
	alpha_1=(RA/(RA+K_r1));

	// <component name="law1_alpha2">
	alpha_2=(F/(F+K_r2));
}