/*
 * Guyton Model: Circulatory Dynamics
 * 
 * Model Status
 * 
 * This CellML model has not been validated. The equations in this
 * file may contain errors and the output from the model may not
 * conform to the results from the MODSIM program. Due to the differences
 * between procedural code (in this case C-code) and declarative
 * languages (CellML), some aspects of the original model were
 * not able to be encapsulated by the CellML model (such as the
 * damping of variables). Work is underway to fix these omissions
 * and validate the CellML model. We also anticipate that many
 * of these problems will be fixed when the CellML 1.0 models are
 * combined in a CellML 1.1 format.
 * 
 * Model Structure
 * 
 * Arthur Guyton (1919-2003) was an American physiologist who became
 * famous for his 1950s experiments in which he studied the physiology
 * of cardiac output and its relationship with the peripheral circulation.
 * The results of these experiments challenged the conventional
 * wisdom that it was the heart itself that controlled cardiac
 * output. Instead Guyton demonstrated that it was the need of
 * the body tissues for oxygen which was the real regulator of
 * cardiac output. The "Guyton Curves" describe the relationship
 * between right atrial pressures and cardiac output, and they
 * form a foundation for understanding the physiology of circulation.
 * 
 * The Guyton model of fluid, electrolyte, and circulatory regulation
 * is an extensive mathematical model of human circulatory physiology,
 * capable of simulating a variety of experimental conditions,
 * and contains a number of linked subsystems relating to circulation
 * and its neuroendocrine control.
 * 
 * This is a CellML translation of the Guyton model of the regulation
 * of the circulatory system. The complete model consists of separate
 * modules each of which characterise a separate physiological
 * subsystems. The Circulation Dynamics is the primary system,
 * to which other modules/blocks are connected. The other modules
 * characterise the dynamics of the kidney, electrolytes and cell
 * water, thirst and drinking, hormone regulation, autonomic regulation,
 * cardiovascular system etc, and these feedback on the central
 * circulation model. The CellML code in these modules is based
 * on the C code from the programme C-MODSIM created by Dr Jean-Pierre
 * Montani.
 * 
 * This particular CellML model describes the flow of blood around
 * the circuit from arteries, to veins, to heart, to lungs, and
 * back to heart again. It also calculates the resistances and
 * the effect of various factors on the resistances. In other words,
 * this section presents the basic hemodynamics of the circulatory
 * system.
 * 
 * model diagram
 * 
 * [[Image file: full_model.png]]
 * 
 * A systems analysis diagram for the full Guyton model describing
 * circulation regulation.
 * 
 * model diagram
 * 
 * [[Image file: circ_dyn.png]]
 * 
 * A schematic diagram of the components and processes described
 * in the current CellML model.
 * 
 * There are several publications referring to the Guyton model.
 * One of these papers is cited below:
 * 
 * Circulation: Overall Regulation, A.C. Guyton, T.G. Coleman,
 * and H.J. Granger, 1972, Annual Review of Physiology , 34, 13-44.
 * PubMed ID: 4334846
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 1 equation(s)
unit conversion off;
unit minute=60 second^1;
unit per_minute=.01666667 second^(-1);
//Warning:  unit mmHg_ renamed from mmHg, as the latter is predefined in JSim with different fundamental units.
unit mmHg_=133.322 kilogram^1*meter^(-1)*second^(-2);
unit per_mmHg=.00750064 kilogram^(-1)*meter^1*second^2;
unit mmHg_L=.133322 kilogram^1*meter^2*second^(-2);
unit per_mmHg2=5.6259564E-5 kilogram^(-2)*meter^2*second^4;
unit mmHg3=2.369766E6 kilogram^3*meter^(-3)*second^(-6);
unit monovalent_mEq=.001 mole^1;
unit monovalent_mEq_per_minute=1.6666667E-5 second^(-1)*mole^1;
unit monovalent_mEq_per_litre=1 meter^(-3)*mole^1;
unit monovalent_mEq_per_litre_per_minute=.01666667 meter^(-3)*second^(-1)*mole^1;
unit litre2_per_monovalent_mEq_per_minute=1.6666667E-5 meter^6*second^(-1)*mole^(-1);
unit L_per_minute=1.6666667E-5 meter^3*second^(-1);
unit mL=1E-6 meter^3;
unit gram_per_L=1 kilogram^1*meter^(-3);
unit L_mmHg_per_gram=133.322 meter^2*second^(-2);
unit L2_mmHg_per_gram2=133.322 kilogram^(-1)*meter^5*second^(-2);
unit mmHg_minute_per_L=7999320 kilogram^1*meter^(-4)*second^(-1);
unit gram_per_minute=1.6666667E-5 kilogram^1*second^(-1);
unit mL_per_L=.001  dimensionless;
unit mL_per_L_per_mmHg=7.5006376E-6 kilogram^(-1)*meter^1*second^2;
unit mL_per_L_per_minute=1.6666667E-5 second^(-1);
unit mL_per_minute_per_mmHg=1.2501063E-10 kilogram^(-1)*meter^4*second^1;
unit L_mL_per_minute_per_mmHg=1.2501063E-13 kilogram^(-1)*meter^7*second^1;
unit L_per_mL=1E3  dimensionless;
unit L_per_mmHg=7.5006376E-6 kilogram^(-1)*meter^4*second^2;
unit mL_per_minute=1.6666667E-8 meter^3*second^(-1);
unit L_per_minute_per_mmHg=1.2501063E-7 kilogram^(-1)*meter^4*second^1;
unit L_per_minute_per_mmHg2=9.376594E-10 kilogram^(-2)*meter^5*second^3;
unit mmHg_per_mL=1.33322E8 kilogram^1*meter^(-4)*second^(-2);

math main {
	realDomain time minute;
	time.min=0;
	extern time.max;
	extern time.delta;
	real ADHMV dimensionless;
	ADHMV=1.0;
	real AMM dimensionless;
	AMM=1.0;
	real ANU dimensionless;
	ANU=0.925271;
	real ANUVN dimensionless;
	ANUVN=1.0;
	real ARM dimensionless;
	ARM=1.16463;
	real ATRRFB dimensionless;
	ATRRFB=1.0;
	real ATRVFB litre;
	ATRVFB=0.0;
	real AU dimensionless;
	AU=1.00022;
	real AUH dimensionless;
	AUH=1.00012;
	real AUM dimensionless;
	AUM=1.00066;
	real AVE dimensionless;
	AVE=1.0;
	real HMD dimensionless;
	HMD=1.0;
	real HPL dimensionless;
	HPL=1.00163;
	real HPR dimensionless;
	HPR=1.00237;
	real MYOGRS dimensionless;
	MYOGRS=1.0;
	real OSA dimensionless;
	OSA=0.97287;
	real PAMK dimensionless;
	PAMK=1.0;
	real PC mmHg_;
	PC=16.9144;
	real RBF L_per_minute;
	RBF=1.22057;
	real VIM dimensionless;
	VIM=1.00076;
	real VP litre;
	VP=3.00449;
	real VRC litre;
	VRC=2.00439;
	real VV6 litre;
	VV6=0.0101913;
	real VV7 litre;
	VV7=0.00366525;
	real VVR litre;
	VVR=2.50967;
	real VVS1(time) litre;
	when(time=time.min) VVS1=3.28246;
	real VAS1(time) litre;
	when(time=time.min) VAS1=0.862514;
	real VLA1(time) litre;
	when(time=time.min) VLA1=0.379883;
	real VPA1(time) litre;
	when(time=time.min) VPA1=0.38131;
	real VRA1(time) litre;
	when(time=time.min) VRA1=0.100043;
	real VBD(time) litre;
	real QVO(time) L_per_minute;
	real QRO(time) L_per_minute;
	real VRA(time) litre;
	real DRA(time) L_per_minute;
	real PRA(time) mmHg_;
	real VRE(time) litre;
	real PRA1(time) mmHg_;
	real HTAUML dimensionless;
	HTAUML=0.4;
	real PPA(time) mmHg_;
	real RVM(time) dimensionless;
	real PP2(time) mmHg_;
	real QLO(time) L_per_minute;
	real QLN(time) L_per_minute;
	real HPEF(time) L_per_minute;
	real QRF L_per_minute;
	QRF=0.15;
	real HSR dimensionless;
	HSR=1;
	real QRN(time) dimensionless;
	real QPO(time) L_per_minute;
	real VPA(time) litre;
	real DPA(time) L_per_minute;
	real VPE(time) litre;
	real RPA(time) mmHg_minute_per_L;
	real PP1T(time) L_per_minute_per_mmHg;
	real PP1(time) L_per_minute_per_mmHg;
	real CPA(time) L_per_minute_per_mmHg;
	real PLA(time) mmHg_;
	real RPV(time) mmHg_minute_per_L;
	real PL1(time) mmHg_;
	real RPT(time) mmHg_minute_per_L;
	real PGL(time) mmHg_;
	real VLA(time) litre;
	real DLA(time) L_per_minute;
	real VLE(time) litre;
	real PLA1(time) mmHg_;
	real PA(time) mmHg_;
	real LVM(time) dimensionless;
	real PA2(time) mmHg_;
	real QLOT(time) L_per_minute;
	real HSL dimensionless;
	HSL=1;
	real QLO1(time) L_per_minute;
	real QAO(time) L_per_minute;
	real VVS(time) litre;
	real DVS(time) L_per_minute;
	real VVA litre;
	real ANY litre;
	ANY=-0.2;
	real VVE(time) litre;
	real VVE1(time) litre;
	real PVS(time) mmHg_;
	real CV L_per_mmHg;
	CV=0.1;
	real PVS1(time) mmHg_;
	real PR1(time) mmHg_;
	real PR1LL mmHg_;
	PR1LL=0;
	real RVG(time) mmHg_minute_per_L;
	real PGV(time) mmHg_;
	real RVS mmHg_minute_per_L;
	real CN7 dimensionless;
	CN7=0.2;
	real CN2 per_mmHg;
	CN2=0.0212;
	real RVSM mmHg_minute_per_L;
	RVSM=1;
	real CN3 dimensionless;
	real RV1 mmHg_minute_per_L;
	real NNRVR mmHg_minute_per_L;
	real VAS(time) litre;
	real DAS(time) L_per_minute;
	real PAG(time) mmHg_;
	real VAE(time) litre;
	real PAM(time) dimensionless;
	real PAEX dimensionless;
	PAEX=2;
	real R1(time) dimensionless;
	real NNRAR(time) mmHg_minute_per_L;
	real RAR mmHg_minute_per_L;
	RAR=30.52;
	real RMULT1 dimensionless;
	RMULT1=1;
	real PGS(time) mmHg_;
	real RSM(time) mmHg_minute_per_L;
	real RAM mmHg_minute_per_L;
	RAM=96.3;
	real RSN(time) mmHg_minute_per_L;
	real BFM(time) L_per_minute;
	real BFN(time) L_per_minute;
	real FISFLO(time) L_per_minute;
	real FIS L_per_minute_per_mmHg;
	FIS=0;
	real SYSFLO(time) L_per_minute;
	real RTP(time) mmHg_minute_per_L;

	// <component name="environment">

	// <component name="circulatory_dynamics">

	// <component name="total_blood_volume_change">
	VBD=((VP+VRC-VVS1-VAS1-VLA1-VPA1-VRA1)/2);

	// <component name="right_atrium">

	// <component name="right_atrial_blood_volume">
	DRA=(QVO-QRO);
	VRA1:time=DRA;
	VRA=(VRA1+VBD*.0574);

	// <component name="right_atrial_pressure">
	VRE=(VRA-(.1 litre));
	PRA=(VRE/(.005 L_per_mmHg));

	// <component name="autonomic_stimulation_effect_on_right_atrial_pressure">
	PRA1=((PRA+(8 mmHg_))*(HTAUML*(AU-1)+1)-(8 mmHg_));

	// <component name="right_ventricle">

	// <component name="pressure_effect_on_right_ventricular_pumping">
	PP2=(PPA/AUH/OSA);
	RVM=(if (PP2<=(0 mmHg_)) 1.06 else if ((PP2>(0 mmHg_)) and (PP2<=(32 mmHg_))) 1.06+((.97 per_mmHg)-(1.06 per_mmHg))*(PP2-(0 mmHg_))/(32-0) else if ((PP2>(32 mmHg_)) and (PP2<=(38.4 mmHg_))) .97+((.93 per_mmHg)-(.97 per_mmHg))*(PP2-(32 mmHg_))/(38.4-32) else if ((PP2>(38.4 mmHg_)) and (PP2<=(48 mmHg_))) .93+((.8 per_mmHg)-(.93 per_mmHg))*(PP2-(38.4 mmHg_))/(48-38.4) else if ((PP2>(48 mmHg_)) and (PP2<=(60.8 mmHg_))) .8+((.46 per_mmHg)-(.8 per_mmHg))*(PP2-(48 mmHg_))/(60.8-48) else if ((PP2>(60.8 mmHg_)) and (PP2<=(72 mmHg_))) .46+((0 per_mmHg)-(.46 per_mmHg))*(PP2-(60.8 mmHg_))/(72-60.8) else 0);

	// <component name="pumping_effectiveness_of_right_ventricle">
	HPEF=(((1 L_per_minute)-QRF)*AUH*RVM*HSR*HMD*HPR+QRF*QLO/QLN);

	// <component name="right_ventricular_output">
	QRN=(if (PRA1<=((-1)*(8 mmHg_))) 0 else if ((PRA1>((-1)*(8 mmHg_))) and (PRA1<=((-1)*(6 mmHg_)))) 0+((.75 per_mmHg)-(0 per_mmHg))*(PRA1-(-1)*(8 mmHg_))/((-1)*6-(-1)*8) else if ((PRA1>((-1)*(6 mmHg_))) and (PRA1<=((-1)*(2 mmHg_)))) .75+((2.6 per_mmHg)-(.75 per_mmHg))*(PRA1-(-1)*(6 mmHg_))/((-1)*2-(-1)*6) else if ((PRA1>((-1)*(2 mmHg_))) and (PRA1<=(4 mmHg_))) 2.6+((9.8 per_mmHg)-(2.6 per_mmHg))*(PRA1-(-1)*(2 mmHg_))/(4-(-1)*2) else if ((PRA1>(4 mmHg_)) and (PRA1<=(12 mmHg_))) 9.8+((13.5 per_mmHg)-(9.8 per_mmHg))*(PRA1-(4 mmHg_))/(12-4) else 13.5);
	QRO=(QRN*HPEF);

	// <component name="pulmonary_vasculature">

	// <component name="pulmonary_vasculature_blood_volume">
	DPA=(QRO-QPO);
	VPA1:time=DPA;
	VPA=(VPA1+VBD*.155);

	// <component name="pulmonary_vasculature_pressure">
	VPE=(VPA-(.30625 litre));
	PPA=(VPE/(.0048 L_per_mmHg));

	// <component name="pulmonary_arterial_resistance">
	PP1T=((.026 L_per_minute_per_mmHg2)*PPA);
	PP1=(if (PP1T<(1E-5 L_per_minute_per_mmHg)) (1E-5 L_per_minute_per_mmHg) else PP1T);
	CPA=(PP1^.5);
	RPA=(1/CPA);

	// <component name="pulmonary_venous_resistance">
	PL1=(PLA+(18 mmHg_));
	RPV=(1/(PL1*(.0357 L_per_minute_per_mmHg2)));

	// <component name="total_pulmonary_vascular_resistance">
	RPT=(RPV+RPA);

	// <component name="pressure_gradient_through_the_lungs">
	PGL=(PPA-PLA);

	// <component name="rate_of_blood_flow_from_pulmonary_veins_to_left_atrium">
	QPO=(PGL/RPT);

	// <component name="left_atrium">

	// <component name="left_atrial_blood_volume">
	DLA=(QPO-QLO);
	VLA1:time=DLA;
	VLA=(VLA1+VBD*.128);

	// <component name="left_atrial_pressure">
	VLE=(VLA-(.38 litre));
	PLA=(VLE/(.01 L_per_mmHg));

	// <component name="autonomic_stimulation_effect_on_left_atrial_pressure">
	PLA1=((PLA+(4 mmHg_))*(HTAUML*(AU-1)+1)-(4 mmHg_));

	// <component name="left_ventricle">

	// <component name="pumping_effectiveness_of_left_ventricle">
	PA2=(PA/(AUH*OSA));
	LVM=(if (PA2<=(0 mmHg_)) 1.04 else if ((PA2>(0 mmHg_)) and (PA2<=(60 mmHg_))) 1.04+((1.025 per_mmHg)-(1.04 per_mmHg))*(PA2-(0 mmHg_))/(60-0) else if ((PA2>(60 mmHg_)) and (PA2<=(125 mmHg_))) 1.025+((.97 per_mmHg)-(1.025 per_mmHg))*(PA2-(60 mmHg_))/(125-60) else if ((PA2>(125 mmHg_)) and (PA2<=(160 mmHg_))) .97+((.88 per_mmHg)-(.97 per_mmHg))*(PA2-(125 mmHg_))/(160-125) else if ((PA2>(160 mmHg_)) and (PA2<=(200 mmHg_))) .88+((.59 per_mmHg)-(.88 per_mmHg))*(PA2-(160 mmHg_))/(200-160) else if ((PA2>(200 mmHg_)) and (PA2<=(240 mmHg_))) .59+((0 per_mmHg)-(.59 per_mmHg))*(PA2-(200 mmHg_))/(240-200) else 0);

	// <component name="left_ventricular_output">
	QLN=(if (PLA1<=((-1)*(2 mmHg_))) (.01 L_per_minute) else if ((PLA1>((-1)*(2 mmHg_))) and (PLA1<=(1 mmHg_))) (.01 L_per_minute)+((3.6 L_per_minute_per_mmHg)-(.01 L_per_minute_per_mmHg))*(PLA1-(-1)*(2 mmHg_))/(1-(-1)*2) else if ((PLA1>(1 mmHg_)) and (PLA1<=(5 mmHg_))) (3.6 L_per_minute)+((9.4 L_per_minute_per_mmHg)-(3.6 L_per_minute_per_mmHg))*(PLA1-(1 mmHg_))/(5-1) else if ((PLA1>(5 mmHg_)) and (PLA1<=(8 mmHg_))) (9.4 L_per_minute)+((11.6 L_per_minute_per_mmHg)-(9.4 L_per_minute_per_mmHg))*(PLA1-(5 mmHg_))/(8-5) else if ((PLA1>(8 mmHg_)) and (PLA1<=(12 mmHg_))) (11.6 L_per_minute)+((13.5 L_per_minute_per_mmHg)-(11.6 L_per_minute_per_mmHg))*(PLA1-(8 mmHg_))/(12-8) else (13.5 L_per_minute));
	QLOT=(LVM*QLN*AUH*HSL*HMD*HPL);
	QLO1=((PLA-PA)/(3 mmHg_minute_per_L));
	QLO=(if (QLO1>(0 L_per_minute)) QLOT+QLO1 else QLOT);

	// <component name="systemic_venous_system">

	// <component name="venous_blood_volume">
	DVS=(QAO-QVO);
	VVS1:time=DVS;
	VVS=(VVS1+VBD*.3986);

	// <component name="angiotensin_induced_venous_constriction">
	VVA=((ANU-1)*ANY);

	// <component name="venous_excess_volume">
	VVE1=(VVS-VVR-VVA-VV7-VV6-ATRVFB);
	VVE=(if (VVE1<(1E-4 litre)) (1E-4 litre) else VVE1);

	// <component name="venous_average_pressure">
	PVS1=((3.7 mmHg_)+(VVE-(.74 litre))/CV);
	PVS=(if (PVS1<(1E-4 mmHg_)) (1E-4 mmHg_) else PVS1);

	// <component name="venous_outflow_pressure_into_heart">
	PR1=(if (PRA<PR1LL) PR1LL else PRA);

	// <component name="resistance_from_veins_to_right_atrium">
	RVG=((.74 mmHg_minute_per_L)/(PVS/(VIM*(3.7 mmHg_)))^.5);

	// <component name="rate_of_blood_flow_from_veins_to_right_atrium">
	PGV=(PVS-PR1);
	QVO=(PGV/RVG);

	// <component name="venous_resistance">
	CN3=(((PC-(17 mmHg_))*CN7+(17 mmHg_))*CN2);
	RV1=(RVSM/CN3);
	RVS=(AVE*RV1*VIM*ANUVN);

	// <component name="NM_NR_venous_resistance">
	NNRVR=(RVS*1.79);

	// <component name="systemic_arterial_system">

	// <component name="arterial_blood_volume">
	DAS=(QLO-QAO);
	VAS1:time=DAS;
	VAS=(VAS1+VBD*.261);

	// <component name="arterial_pressure_and_pressure_gradient">
	VAE=(VAS-(.495 litre));
	PA=(VAE/(.00355 L_per_mmHg));
	PAG=(PA-PRA);

	// <component name="pressure_effect_on_arterial_distention">
	PAM=((PA/(100 mmHg_))^PAEX);

	// <component name="non_renal_systemic_arterial_resistance_multiplier">
	R1=(ANU*ADHMV*AUM*VIM*PAMK/PAM/ATRRFB);

	// <component name="NM_NR_arterial_resistance">
	NNRAR=(RAR*ARM*R1*MYOGRS*RMULT1);

	// <component name="pressure_gradient_from_arteries_to_veins">
	PGS=(PA-PVS);

	// <component name="M_systemic_resistance">
	RSM=(RAM*AMM*R1*MYOGRS*RMULT1);

	// <component name="total_NM_NR_systemic_resistance">
	RSN=(NNRAR+NNRVR);

	// <component name="blood_flow_through_M_tissues">
	BFM=(PGS/RSM);

	// <component name="blood_flow_through_NM_NR_tissues">
	BFN=(PGS/RSN);

	// <component name="blood_flow_through_AV_fistulas">
	FISFLO=(PAG*FIS);

	// <component name="systemic_blood_flow">
	SYSFLO=(BFM+BFN+RBF);
	QAO=(SYSFLO+FISFLO);

	// <component name="total_peripheral_resistance">
	RTP=(PAG/QAO);

	// <component name="parameter_values">
}