tests/testthat/test-issue-117.R

In nlmixrdevelopment/RxODE: Facilities for Simulating from ODE-Based Models

rxodeTest(
  {
    context("Test large model compiles")

    mod <- RxODE("
##### define heart failure parameters. Allow them to increase from baseline over time #####
if (Heart_failure_link == 1) {
R_art0		  = (R_art0_initial-R_art0_initial*heart_failure_resistance_aorta_scale)/exp(sim_time/(24*5))+R_art0_initial*heart_failure_resistance_aorta_scale;
R_per0	  	  = (R_per0_initial-R_per0_initial*heart_failure_resistance_prepheral_scale)/exp(sim_time/(24*5))+R_per0_initial*heart_failure_resistance_prepheral_scale;
R_al0 = R_al0_scale*R_per0;    #0.66*R_per0
R_cap0= R_cap0_scale*R_per0;    #0.25*R_per0
R_vn0 = R_vn0_scale*R_per0;   #0.09*R_per0
R_ven0	  	= (R_ven0_initial-R_ven0_initial*heart_failure_resistance_venous_scale)/exp(sim_time/(24*5))+R_ven0_initial*heart_failure_resistance_venous_scale;
R_art_pulm	= (R_art_pulm_initial-R_art_pulm_initial*heart_failure_resistance_R_art_pulm_scale)/exp(sim_time/(24*5))+R_art_pulm_initial*heart_failure_resistance_R_art_pulm_scale;
R_ven_pulm	= (R_ven_pulm_initial-R_ven_pulm_initial*heart_failure_resistance_R_ven_pulm_scale)/exp(sim_time/(24*5))+R_ven_pulm_initial*heart_failure_resistance_R_ven_pulm_scale;

C_art	    	= (C_art_initial-C_art_initial*heart_failure_compliance_aorta_scale)/exp(sim_time/(24*5))+C_art_initial*heart_failure_compliance_aorta_scale;
C_ven0	  	= (C_ven0_initial-C_ven0_initial*heart_failure_compliance_venous_scale)/exp(sim_time/(24*5))+C_ven0_initial*heart_failure_compliance_venous_scale;
C_pulm_ven  = (C_pulm_ven_initial-C_pulm_ven_initial*heart_failure_compliance_pulm_ven_scale)/exp(sim_time/(24*5))+C_pulm_ven_initial*heart_failure_compliance_pulm_ven_scale;
C_pulm_art  = (C_pulm_art_initial-C_pulm_art_initial*heart_failure_compliance_pulm_art_scale)/exp(sim_time/(24*5))+C_pulm_art_initial*heart_failure_compliance_pulm_art_scale;

cf            = (cf_initial-cf_initial*cf_scale)/exp(sim_time/(24))+cf_initial*cf_scale;
contractility = (contractility_initial-contractility_initial*contractility_HF_scale)/exp(sim_time/(24))+contractility_initial*contractility_HF_scale;

Sodium_protein_filtration_rate_Kf = (Sodium_protein_filtration_rate_Kf_nom-Sodium_protein_filtration_rate_Kf_nom*kf_scale)/exp(sim_time/(24))+Sodium_protein_filtration_rate_Kf_nom*kf_scale;
vascular_responsiveness_scale =  (vascular_responsiveness_scale_nom-vascular_responsiveness_scale_nom*vascular_responsiveness_scale_nom_scale)/exp(sim_time/(24))+vascular_responsiveness_scale_nom*vascular_responsiveness_scale_nom_scale;
}
else {  #If no heart failure, parameters just equation baseline values
R_art0		  = R_art0_initial	;
R_per0		  = R_per0_initial	;
R_al0 = R_al0_scale*R_per0;    #0.66*R_per0
R_cap0= R_cap0_scale*R_per0;    #0.25*R_per0
R_vn0 = R_vn0_scale*R_per0;   #0.09*R_per0
R_ven0		  = R_ven0_initial	;
R_art_pulm	= R_art_pulm_initial;
R_ven_pulm	= R_ven_pulm_initial;
C_art		    = C_art_initial 	;
C_ven0	  	= C_ven0_initial	;
C_pulm_ven  = C_pulm_ven_initial;
C_pulm_art  = C_pulm_art_initial;

cf = cf_initial;
contractility = contractility_initial;
Sodium_protein_filtration_rate_Kf = Sodium_protein_filtration_rate_Kf_nom;
vascular_responsiveness_scale = vascular_responsiveness_scale_nom;
}
CO_nom = CO_nom_initial;

plasma_protein_concentration =plasma_protein_amount / (blood_volume_L * L_dL);
ISF_protein_concentration = ISF_protein_amount / (interstitial_fluid_volume * L_dL);


if(test_ISF_pressure==1){
ISF_pressure =  -0.000000002316554*interstitial_fluid_volume^6
+ 0.000000910096686*interstitial_fluid_volume^5 - 0.000141950241878*interstitial_fluid_volume^4
+ 0.011247538513374*interstitial_fluid_volume^3 - 0.471694546620885*interstitial_fluid_volume^2
+ 9.92251133652023*interstitial_fluid_volume - 81.2426416630404+5.09732;

}else{
ISF_pressure=ISF_pressure_initial;
}

########### feedback of the cardiac output perturbation#################

#If not testing cardiac reserve:
if (Heart_failure_link == 0 & cardiac_reserve == 0  ) {
tissue_autoregulation_signal = max(0.1,1+tissue_autoreg_scale*(Kp_CO*(CO_delayed - CO_nom*CO_species_scale)+Ki_CO*CO_error));
peripheral_resistance_multiplier = disease_effect_on_TPR_peripheral_resistance * B2sna_effect_on_TPR*A1sna_effect_on_TPR*tissue_autoregulation_signal;
peripheral_resistance_multiplier_adjusted = max(dilation_scale,time_TPR_scale*(1+vascular_responsiveness_scale*(peripheral_resistance_multiplier-1)));#*max(0,(-peripheral_resistance_multiplier_adjusted+dilation_scale)/dilation_scale);
heart_rate = HR_heart_rate * BB_HR_effect ;

} else if (Heart_failure_link == 0 & cardiac_reserve == 1){ #If testing cardiac reserve

tissue_autoregulation_signal = max(0.1,1+tissue_autoreg_scale*(Kp_CO*(CO_delayed - CO_nom*CO_species_scale)+Ki_CO*CO_error));
HR_autoregulation_signal = 1+tissue_autoreg_scale*(HRP_CO*(CO_delayed - CO_nom*CO_species_scale)+HRI_CO*CO_error);

heart_rate_multiplier_adjusted = time_HR_scale*(1-HR_autoregulation_signal);
Heart_rate_increasing_ratio =  min(Mag_HR_changing_ratio,(Heart_rate_increasing_ratio + heart_rate_multiplier_adjusted));#*max(0,(-heart_rate_multiplier_adjusted+Mag_HR_changing_ratio)/Mag_HR_changing_ratio);
heart_rate = (HR_heart_rate * BB_HR_effect*(1+Heart_rate_increasing_ratio));

peripheral_resistance_multiplier = disease_effect_on_TPR_peripheral_resistance * B2sna_effect_on_TPR*A1sna_effect_on_TPR*tissue_autoregulation_signal;
peripheral_resistance_multiplier_adjusted = max(dilation_scale,time_TPR_scale*(1+vascular_responsiveness_scale*(peripheral_resistance_multiplier-1)));#*max(0,(-peripheral_resistance_multiplier_adjusted+dilation_scale)/dilation_scale);

} else if (Heart_failure_link == 1 & cardiac_reserve == 0){

tissue_autoregulation_signal = max(0.1,1+tissue_autoreg_scale*(Kp_HF_CO*(CO_delayed - CO_nom*CO_species_scale)+Ki_HF_CO*CO_error));
HR_autoregulation_signal = 1+tissue_autoreg_scale*(HRP_HF_CO*(CO_delayed - CO_nom*CO_species_scale)+HRI_HF_CO*CO_error);

peripheral_resistance_multiplier = disease_effect_on_TPR_peripheral_resistance * B2sna_effect_on_TPR*A1sna_effect_on_TPR*tissue_autoregulation_signal;
peripheral_resistance_multiplier_adjusted =max(dilation_HF_Chronic_scale,(1+vascular_responsiveness_scale*(peripheral_resistance_multiplier-1)));#*max(0,(-peripheral_resistance_multiplier_adjusted+dilation_HF_scale)/dilation_HF_scale);

heart_rate_multiplier_adjusted = (1-HR_autoregulation_signal);
Heart_rate_increasing_ratio = min(Mag_HR_HF_Chronic_ratio,(Heart_rate_increasing_ratio + heart_rate_multiplier_adjusted));#*max(0,(-heart_rate_multiplier_adjusted+Mag_HR_HF_changing_ratio)/Mag_HR_HF_changing_ratio);
heart_rate = (HR_heart_rate * BB_HR_effect*(1+Heart_rate_increasing_ratio));
} else if (Heart_failure_link == 1& cardiac_reserve == 1){

tissue_autoregulation_signal = max(0.1,1+tissue_autoreg_scale*(Kp_HF_CO*(CO_delayed - CO_nom*CO_species_scale)+Ki_HF_CO*CO_error));
HR_autoregulation_signal = 1+tissue_autoreg_scale*(HRP_HF_CO*(CO_delayed - CO_nom*CO_species_scale)+HRI_HF_CO*CO_error);

peripheral_resistance_multiplier = disease_effect_on_TPR_peripheral_resistance * B2sna_effect_on_TPR*A1sna_effect_on_TPR*tissue_autoregulation_signal;
peripheral_resistance_multiplier_adjusted =max(dilation_HF_Acute_scale,(1+vascular_responsiveness_scale*(peripheral_resistance_multiplier-1)));#*max(0,(-peripheral_resistance_multiplier_adjusted+dilation_HF_scale)/dilation_HF_scale);

heart_rate_multiplier_adjusted = (1-HR_autoregulation_signal);
Heart_rate_increasing_ratio = min(Mag_HR_HF_Acute_ratio,(Heart_rate_increasing_ratio + heart_rate_multiplier_adjusted));#*max(0,(-heart_rate_multiplier_adjusted+Mag_HR_HF_changing_ratio)/Mag_HR_HF_changing_ratio);
heart_rate = (HR_heart_rate * BB_HR_effect*(1+Heart_rate_increasing_ratio));
}
########### feedback of the cardiac output perturbation #################


peripheral_resistance =TPR_scale_peripheral_resistance * R_al0*peripheral_resistance_multiplier_adjusted;

arterial_dis_resistance=peripheral_resistance;
capillary_resistance=R_cap0;
venules_resistance=R_vn0;


beat_duration = min_sec / heart_rate ;
beat_time = sim_time/beat_duration - floor(sim_time/beat_duration);
periods = floor(sim_time/beat_duration);

blood_volume = blood_volume_L/1000;
time_step = 0.01;

#Neurohormonal effects - right now they are all turned off   SNA''sympathetic nerve activity''
sna_effect_on_contractility=1;
sna_effect_on_HR = 1;
AngII_effect_on_venous_compliance=1;


# ANP_effect_on_venous_compliance=1;
SNA_effect_on_venous_compliance=1;
B2sna_effect_on_TPR = 1;
A1sna_effect_on_TPR =  1;
Ang_II_effect_on_systemic_resistance =1;
aldo_effect_on_systemic_resistances = 1;
CCB_effect_on_systemic_arterial_resistance = 1;


# ANP_effect_on_peripheral_resistance = 1;
glu_eff_1 = 0;
angII_eff_1 = 0;
aldo_eff_1 = 0;

#Beta-blockers are drugs that bind to beta-adrenoceptors
#and thereby block the binding of norepinephrine and epinephrine to these receptors.
beta_blocker_effect_on_contractility = BB_contractility_effect; # BB_contractility_effect = 1
BP_effect_on_compliance=1;




##### Volume unit conversions #####
LV_volume_mL = LV_volume * m3_mL;
arterial_volume_mL = arterial_volume * m3_mL;
arterial_dis_circulation_volume_mL 	=arterial_dis_circulation_volume *m3_mL;
capillary_circulation_volume_mL	 	=capillary_circulation_volume	 *m3_mL;
venules_circulation_volume_mL		=venules_circulation_volume		 *m3_mL;

############

RV_volume_mL = RV_volume * m3_mL;
pulmonary_arterial_volume_mL = pulmonary_arterial_volume * m3_mL;
venous_volume_mL = venous_volume * m3_mL;
#total_blood_volume_mL = LV_volume_mL + arterial_volume_mL + peripheral_volume_mL + RV_volume_mL + pulmonary_arterial_volume_mL + pulmonary_venous_volume * m3_mL + venous_volume_mL;
total_blood_volume_mL = LV_volume_mL + arterial_volume_mL + arterial_dis_circulation_volume_mL+capillary_circulation_volume_mL+venules_circulation_volume_mL + RV_volume_mL + pulmonary_arterial_volume_mL + pulmonary_venous_volume * m3_mL + venous_volume_mL;



########## Cardiac tissue composition ################
# V_w_0 = 0.00012
# Baseline_Interstitial_Fibrosis = 0.000003
#Baseline_Replacement_Fibrosis = 0.000003
#Baseline_Interstitial_Tissue = 0.000032
#Baseline_Segmental_Fibrosis = 0

baseline_total_myocyte_volume = V_w_0 - Baseline_Interstitial_Fibrosis - Baseline_Replacement_Fibrosis - Baseline_Interstitial_Tissue;  ## baseline myocyte volume determined by V_w_0
baseline_single_myocyte_volume = baseline_total_myocyte_volume/Baseline_Myocyte_Number;
baseline_myocyte_diameter = 2*sqrt(baseline_single_myocyte_volume/(Pi*Baseline_Myocyte_Length));	## baseline myocyte diameter determined by V_w_0 & Baseline_Myocyte_Length, NOT by Baseline_Myocyte_Diameter

myocyte_length = Baseline_Myocyte_Length + change_in_myocyte_length;						## change_in_myocyte_length depends on passive stress levels
myocyte_diameter = baseline_myocyte_diameter + change_in_myocyte_diameter;					## change_in_myocyte_diameter depends on active stress levels - HYPERTROPHY
single_myocyte_volume = myocyte_length * Pi * (myocyte_diameter ^ 2) / 4;					## myocyte volume calculated as a cylinder

#Baseline_Myocyte_Number = 3.3e+9
number_of_live_myocytes = Baseline_Myocyte_Number;
total_myocyte_volume = single_myocyte_volume * number_of_live_myocytes;
total_nonmyocyte_volume = Baseline_Interstitial_Fibrosis + Baseline_Interstitial_Tissue + Baseline_Replacement_Fibrosis;
LV_wall_volume = total_myocyte_volume + total_nonmyocyte_volume;


level_of_hypertrophy = LV_wall_volume / (baseline_total_myocyte_volume + total_nonmyocyte_volume);		## a measure of how much LV wall wolume has grown
outward_growth = LV_cavity_volume / LV_V0_baseline;										## a measure of how much the LV chamber has grown in volume

pct_change_in_myocyte_diameter = change_in_myocyte_diameter / baseline_myocyte_diameter * 100;
pct_change_in_myocyte_length = change_in_myocyte_length / Baseline_Myocyte_Length * 100;


#################################################################################################
#################################################################################################

######### Cardiac Mechanics #########
## Muscle fiber stress and strain are approximately homogeneously distributed, so that they may be approximated by single values.
## Microscopic constitutive laws for fiber stress and radial stress are used to model active and passive fiber stress.

LV_cavity_volume = LV_V0_baseline * (1 + myo_L_scale * change_in_myocyte_length / Baseline_Myocyte_Length) ^ 3 * (1 - myo_D_scale * change_in_myocyte_diameter / baseline_myocyte_diameter) ^ 2;
LV_fiber_stretch =((LV_volume + (LV_wall_volume/3)) / (LV_cavity_volume + (LV_wall_volume / 3)))^0.3333333;

LV_sarcomere_length = ls_0_passive_LV_sarcomere_length * LV_fiber_stretch;
LV_sarcomere_contraction_velocity = (LV_sarcomere_length - LV_sarcomere_length_delayed) / time_step;
contraction_velocity_effect_in_LV = (1 - LV_sarcomere_contraction_velocity / v0_LV_contraction_velocity_effect_in_LV) / (1 + Cv_contraction_velocity_effect_in_LV * LV_sarcomere_contraction_velocity / v0_LV_contraction_velocity_effect_in_LV);

if (LV_sarcomere_length > ls_a0) {
sarcomere_length_effect_in_LV = ((LV_sarcomere_length - ls_a0) / (ls_ar_sarcomere_length_effect_in_LV - ls_a0));
} else {
sarcomere_length_effect_in_LV = 0;
}

chamber_radius = ((LV_cavity_volume * 3 / 4 / Pi) ^ 0.3333333) * m_mm ;			## approximating the LV as a spherical shell
chamber_diameter = 2 * chamber_radius;

outer_radius = (((LV_cavity_volume + LV_wall_volume) * 3 / 4 / Pi) ^ 0.3333333) * m_mm ;
h_wall = outer_radius - chamber_radius;			## wall thickness

h_over_r = h_wall / chamber_radius;				## a measure of the LV chamber growth; it's the ratio of wall thickness to chamber radius


EDV_chamber_radius = ((LV_EDV* 3 / 4 / Pi) ^ 0.33333333) * m_mm;
EDV_chamber_diameter = 2*EDV_chamber_radius;

EDV_outer_radius = (((LV_EDV + LV_wall_volume) * 3 / 4 / Pi) ^ 0.3333333) * m_mm ;
EDV_h_wall = EDV_outer_radius - EDV_chamber_radius;			## wall thickness
EDV_h_over_r = EDV_h_wall / EDV_chamber_radius;

LV_mass = 1000000*LV_wall_volume*1.05 ;			##wall volume*[(cm->m)^3]*density
## density= 1.05 g/cm^3
#LV_mass = ( 0.8*(1.04*((chamber_diameter + h_wall + h_wall)^3 - chamber_diameter^3)) + 0.6 ) / 1000;	# LV mass in g

#### Cardiac excitation ####

RV_twitch_duration = RV_systolic_time_fraction * beat_duration;

t_d = tau_d_LV_twitch_shape;
t_r = tau_r_LV_twitch_shape;
t_twitch = t_r + t_d;

if (beat_time <= t_r) {
sin_signal=(sin(Pi*beat_time/t_twitch))^n_r_LV_twitch_shape;
}else {
sin_signal=(sin(Pi*beat_time/t_twitch))^n_r_LV_twitch_shape;
}

LV_twitch_shape = sin_signal;
if (beat_time < 0) {
LV_twitch_shape = 0;
}
if (beat_time > t_twitch) {
LV_twitch_shape=0;
}

RV_twitch_shape = (sin(Pi * beat_time / RV_twitch_duration)) ^ 2;
if (beat_time < 0) {
RV_twitch_shape=0;
}
if (beat_time > RV_twitch_duration) {
RV_twitch_shape = 0;
}

LV_active_stress = contractility_scale_LV_active_stress * contractility * sigma_ar * sarcomere_length_effect_in_LV * LV_twitch_shape * contraction_velocity_effect_in_LV * beta_blocker_effect_on_contractility *sna_effect_on_contractility;

#Hypertrophy increases ventricular stiffness
hypertrophy_effect_on_Cf = hypertrophy_Cf_slope*(level_of_hypertrophy - 1);
C_f = cf*(1+hypertrophy_effect_on_Cf);


stretch_zero_S = stretch_min_LV_passive_stress_along_fiber - stretch_scale_LV_passive_stress_along_fiber;
if (LV_fiber_stretch >= stretch_zero_S) {
LV_passive_stress_along_fiber = s_f0 * (exp(C_f * (LV_fiber_stretch - stretch_zero_S)) - 1);
} else {
LV_passive_stress_along_fiber = 0;
}


LV_radial_stretch = 1/ (LV_fiber_stretch * LV_fiber_stretch);

if (LV_radial_stretch >= 1) {
LV_passive_radial_stress = s_r0 * (exp(c_r_LV * (LV_radial_stretch - 1)) - 1);
} else {
LV_passive_radial_stress = 0;
}

LV_total_stress = LV_active_stress + LV_passive_stress_along_fiber - 2 * LV_passive_radial_stress;


if (LV_volume > LV_V0_min) {
rel_volume_LV = 1+LV_wall_volume/LV_volume;
} else {
rel_volume_LV = 1+LV_wall_volume/LV_V0_min;
}

LV_pressure = LV_total_stress * log(rel_volume_LV)/ 3;



#################################################################################################
#################################################################################################

#peripheral_pressure = P_ven0 + (peripheral_circulation_volume - V_per0) / C_per;
arterial_dis_pressure = P_al0 + (arterial_dis_circulation_volume-V_al0)/C_al;
capillary_pressure = P_cap0 + (capillary_circulation_volume-V_cap0)/C_cap;
venules_pressure = P_vn0 + (venules_circulation_volume-V_vn0)/C_vn;

# venous_compliance=C_ven0*AngII_effect_on_venous_compliance*ANP_effect_on_venous_compliance*SNA_effect_on_venous_compliance;
venous_compliance=C_ven0*AngII_effect_on_venous_compliance*SNA_effect_on_venous_compliance;
venous_pressure = P_ven0 + (venous_volume - V_ven0) / venous_compliance; ################################################################################# was P_ven0

RV_Cavity_Volume = RV_V0;
RV_wall_volume = V_w_0_RV;
RV_fiber_stretch = ((RV_volume + V_w_0_RV/3) / (RV_Cavity_Volume + RV_wall_volume/3))^(0.333);
RV_sarcomere_length = ls_a0_RV * RV_fiber_stretch;

if (RV_sarcomere_length > ls_a0_RV) {
sarcomere_length_effect_in_RV = (RV_sarcomere_length - ls_a0_RV) / (0.000002 - ls_a0_RV);
} else {
sarcomere_length_effect_in_RV = 0;
}


RV_sarcomere_contraction_velocity = (RV_sarcomere_length - RV_sarcomere_length_delayed) / time_step;
contraction_velocity_effect_in_RV = (1 - RV_sarcomere_contraction_velocity / v0_RV_contraction_velocity_effect_in_RV) / (1 + 0 * RV_sarcomere_contraction_velocity / v0_RV_contraction_velocity_effect_in_RV);
RV_active_stress_multiplier = contractility_RV*sna_effect_on_contractility;
RV_active_stress = contractility_RV * sigma_ar_RV * sarcomere_length_effect_in_RV * RV_twitch_shape * contraction_velocity_effect_in_RV*sna_effect_on_contractility;

RV_radial_stretch = 1/ (RV_fiber_stretch * RV_fiber_stretch);

if (RV_radial_stretch >= 1) {
RV_passive_radial_stress = s_r0_RV * (exp(c_r_RV * (RV_radial_stretch - 1)) - 1);
} else {
RV_passive_radial_stress = 0;
}


if (RV_fiber_stretch >= 1) {
RV_passive_stress_along_fiber = s_f0_RV * (exp(cf_RV * (RV_fiber_stretch - 1)) - 1);
} else {
RV_passive_stress_along_fiber = 0;
}

RV_total_stress = RV_active_stress + RV_passive_stress_along_fiber - 2 * RV_passive_radial_stress;

if (RV_volume > RV_V0_min) {
rel_volume = (1 + RV_wall_volume / RV_volume);
} else {
rel_volume = (1 + RV_wall_volume / RV_V0_min);
}

RV_pressure = RV_total_stress * log(rel_volume) / 3;

######### Circulatory hemodynamics #########
venous_flow = (venules_pressure - venous_pressure) / (R_ven0);

##allow blood volume link between heart and kidney to be turned on/off
if (heart_renal_link == 1) {
venous_volume_target =  blood_volume - LV_volume - arterial_volume - arterial_dis_circulation_volume-capillary_circulation_volume-venules_circulation_volume - RV_volume - pulmonary_arterial_volume - pulmonary_venous_volume;
} else {
venous_volume_target = venous_volume;
}

tricuspid_valve_flow_rate = max((venous_pressure - RV_pressure) / R_r_atrium,min_flux);

pulmonary_arterial_pressure = ( pulmonary_arterial_volume - V_pulm_art0 )/C_pulm_art + P_art0;
pulmonary_venous_pressure = P_ven0 + (pulmonary_venous_volume  - V_pulm_ven0 )/C_pulm_ven;
pulmonary_arterial_blood_flow = (pulmonary_arterial_pressure  - pulmonary_venous_pressure )/ R_ven_pulm ;

dP = RV_pressure - pulmonary_arterial_pressure;
Zn = L_pulm + time_step * R_art_pulm;
pulmonary_blood_flow = (pulmonary_blood_flow_delayed * L_pulm + dP * time_step) / Zn;

pulmonary_valve_flow_rate = max(pulmonary_blood_flow,min_flux);


#Allow regurgitation if pressure differential across mitral valve exceeds threshold (set threshold really high to eliminate any regurgitation)
if (pulmonary_venous_pressure > LV_pressure) {
mitral_valve_flow_rate = max  (  (pulmonary_venous_pressure - LV_pressure)/R_left_atrium  , min_flux ) ;
} else {
if (LV_pressure - pulmonary_venous_pressure < mitral_regurgitation_pressure_diff) {
mitral_valve_flow_rate = min_flux;
} else {
#Allow negative flow back out of the ventricle
mitral_valve_flow_rate = (pulmonary_venous_pressure - LV_pressure)/R_left_atrium;
}

}

pulmonary_valve_flow_rate = max(pulmonary_blood_flow,min_flux);


## arterial_pressure computed without taking account compliance
## arterial_compliance = compliance_scale_arterial_compliance * (C_art ) * BP_effect_on_compliance;
## arterial_pressure = (arterial_volume - V_art0) / arterial_compliance + P_art0 ;

## arterial_pressure computed taking account compliance
## Constants taken from Safar, M. E., et al. Stiffness of carotid artery wall material and blood pressure in humans: application to antihypertensive therapy and stroke prevention. Stroke 31.3 (2000): 782-790.
## We assume a linear effect as follows.
## BP_effect_on_stiffness = (arterial_pressure - 85)*Stiffness_BP_slope;
## Solving the following system

Stiffness0=1/C_art;
arterial_stiffness = Stiffness0*(1+ (MAP_delayed - nominal_map_setpoint)*Stiffness_BP_slope);
arterial_compliance = 1/arterial_stiffness;
arterial_pressure = (arterial_volume - V_art0) / arterial_compliance + P_art0 ;


#peripheral_pressure = P_ven0 + (peripheral_circulation_volume - V_per0) / C_per;
systemic_blood_flow = (arterial_pressure - arterial_dis_pressure) / arterial_dis_resistance;
arterial_dis_blood_flow = (arterial_dis_pressure-capillary_pressure)/capillary_resistance;
capillary_blood_flow = (capillary_pressure-venules_pressure)/venules_resistance;
venules_blood_flow = (venules_pressure-venous_pressure)/R_ven0;
dP_1 = LV_pressure - arterial_pressure;
Zn_1 = L_art + R_art0 * time_step;
aortic_blood_flow = (aortic_blood_flow_delayed * L_art*L_scale + dP_1 * time_step) / Zn_1;

aortic_valve_flow_rate = max(aortic_blood_flow,min_flux);



#############################################
### Relating BNP and NTP to End Diastolic Stress###

BNP = exp(BNP_factor*((LV_EDS+1736)/5.094)+3.14);
NTP = exp((log(BNP)+1.4919)/1.0694);

###Assume ANP proportional to BNP
# ANP = BNP;

## Relate ANP with other variables following Karaaslan et al. Long-term mathematical model involving renal sympathetic nerve activity, arterial pressure, and sodium excretion. Annals of biomedical engineering 33.11 (2005): 1607-1630.
if (heart_renal_link  == 1) {
# Right atrial pressure (Pra) - Frank-Starling law
Pra = 0.2787*exp(CO_delayed*0.2281); # [cardiac output]=l/min, [Pra]=mmHg
# Normalized atrial natriuretic peptide concentration
Canp = max(0,7.427 - 6.554 / ( 1 + exp( Pra - 3.762) ) + deltaCanp); # For human beings the actual value at normal conditions is about 36 ng/l
lambdaANP = -0.1*Canp+1.1199;
} else {
lambdaANP = 1;
}


#############################################
### Find Diastolic and Systolic values ###

if (beat_time >= (1 - .01/beat_duration) && beat_time < 1) {
LV_pressure_diastolic_max = LV_pressure;
LV_stress_diastolic_max = LV_passive_stress_along_fiber;
LV_volume_maximum = LV_volume;

} else {
LV_pressure_diastolic_max = LV_EDP;
LV_stress_diastolic_max = LV_EDS;
LV_volume_maximum = LV_EDV;
}

LV_EDP_old = LV_pressure_diastolic_max;
LV_EDS_old = LV_stress_diastolic_max;
LV_EDV_old = LV_volume_maximum;



## Method to find Systolic and Diastolic blood pressure
## use current and last two time steps to find local maxima and minima
if (arterial_pressure_delayed < arterial_pressure_bigger_delay) {
systemic_pressure_minimum_1 = arterial_pressure_delayed;
} else {
systemic_pressure_minimum_1 = diastolic_pressure;
}
if (arterial_pressure_delayed < arterial_pressure) {
systemic_pressure_minimum = systemic_pressure_minimum_1;
} else {
systemic_pressure_minimum = diastolic_pressure;
}

if (arterial_pressure_delayed > arterial_pressure_bigger_delay) {
systemic_pressure_maximum_1 = arterial_pressure_delayed;
} else {
systemic_pressure_maximum_1 = systolic_pressure;
}
if (arterial_pressure_delayed > arterial_pressure) {
systemic_pressure_maximum= systemic_pressure_maximum_1;
} else {
systemic_pressure_maximum= systolic_pressure;
}

systolic_pressure_old = systemic_pressure_maximum;
diastolic_pressure_old = systemic_pressure_minimum;



#############################################################
### Find Peak systolic stress in the LV ###

if (beat_time >= (t_r*0.8) && beat_time < (t_r*0.85)) {
LV_peak_stress = LV_active_stress;
} else {
LV_peak_stress = LV_active_stress_peak;
}

## fixes a problem where it was finding peaks that were very close to zero
if (LV_active_stress > 1) {
LV_active_stress_peak_old = LV_peak_stress;
} else {
LV_active_stress_peak_old = LV_active_stress_peak;
}


################################Find Capillary and Venous Pressures

###Method to compute the systolic and diastolic venous_pressure
if (venous_pressure_delayed < venous_pressure_bigger_delay) {
venous_pressure_minimum_1 = venous_pressure_delayed;
} else {
venous_pressure_minimum_1 = venous_diastolic_pressure;
}
if (venous_pressure_delayed < venous_pressure) {
venous_pressure_minimum = venous_pressure_minimum_1;
} else {
venous_pressure_minimum = venous_diastolic_pressure;
}

if (venous_pressure_delayed > venous_pressure_bigger_delay) {
venous_pressure_maximum_1 = venous_pressure_delayed;
} else {
venous_pressure_maximum_1 = venous_systolic_pressure;
}
if (venous_pressure_delayed > venous_pressure) {
venous_pressure_maximum= venous_pressure_maximum_1;
} else {
venous_pressure_maximum= venous_systolic_pressure;
}

venous_systolic_pressure_old = venous_pressure_maximum;
venous_diastolic_pressure_old = venous_pressure_minimum;

###Method to compute the systolic and diastolic peripheral_pressure

if (capillary_pressure_delayed < capillary_pressure_bigger_delay) {
capillary_pressure_minimum_1 = capillary_pressure_delayed;
} else {
capillary_pressure_minimum_1 = capillary_diastolic_pressure;
}
if (capillary_pressure_delayed < capillary_pressure) {
capillary_pressure_minimum = capillary_pressure_minimum_1;
} else {
capillary_pressure_minimum = capillary_diastolic_pressure;
}

if (capillary_pressure_delayed > capillary_pressure_bigger_delay) {
capillary_pressure_maximum_1 = capillary_pressure_delayed;
} else {
capillary_pressure_maximum_1 = capillary_systolic_pressure;
}
if (capillary_pressure_delayed > capillary_pressure) {
capillary_pressure_maximum= capillary_pressure_maximum_1;
} else {
capillary_pressure_maximum= capillary_systolic_pressure;
}

capillary_systolic_pressure_old = capillary_pressure_maximum;
capillary_diastolic_pressure_old = capillary_pressure_minimum;



#Calculate mean pressures

mean_capillary_pressure = (capillary_systolic_pressure/3+capillary_diastolic_pressure*2/3)*Pa_mmHg;



if (heart_renal_link  == 1) {
mean_arterial_pressure_MAP = (systolic_pressure/3+diastolic_pressure*2/3)*Pa_mmHg;
mean_venous_pressure = (venous_systolic_pressure+venous_diastolic_pressure*2)/3*Pa_mmHg;
} else {
mean_arterial_pressure_MAP = nominal_map_setpoint;
mean_venous_pressure = P_venous;
}




#################################################################################################
#################################################################################################

######################### Hypertrophy ########################################
## growth occurs when the peak of the active stress exceeds the threshhold
## use the same method as above to locate the peaks of the active stress

######################### Hypertrophy ########################################

if (LV_active_stress_peak > LV_active_stress_threshhold) {		## myocyte diameter grows if peak active stress exceeds the threshhold
kD_hypertrophy =  (kD_HYPERTROPHY*C_renal_CV_timescale) * max(0, ((max_myocyte_diameter_increase) - change_in_myocyte_diameter)/(max_myocyte_diameter_increase));	## progression of hypertophic growth - myocyte diameter increase
} else {																							## there is a limit on how much the diameter of the myocyte can increase (right now: 40%)
kD_hypertrophy = (kD_HYPERTROPHY*C_renal_CV_timescale);		## regression of hypertrophic growth (negative will come from ODE below)
}

if (LV_EDS > LV_passive_stress_along_fiber_threshhold) {
kL_hypertrophy =  (kL_HYPERTROPHY*C_renal_CV_timescale) * max(0, ((max_myocyte_length_increase) - change_in_myocyte_length)/(max_myocyte_length_increase));		## progression of growth - myocyte length increase
} else {              ##adding a length threshhold to limit volumetric remodeling
kL_hypertrophy = 0;		## myocyte length does not decrease once stretched
}

#################################################################################################
#################################################################################################

LVID=((6*LV_EDV)/3.14159)^(1/3);

#############################################################################


############################ Kidney Model ###############################################

#Disease effects on nephrons

#For now we limit glomeruli loss to 30% of Kf damage - somewhat arbitrary
number_of_functional_glomeruli = baseline_nephrons*(1 - 0.3*disease_effects_decreasing_Kf);

#Assume that other nephrons are lost due to tubular effects, which do not affect the glomerulus
number_of_functional_tubules = baseline_nephrons*(1-disease_effect_on_nephrons);


###AT1-bound AngII constricts the preafferent, afferent, and efferent arterioles
AT1_preaff_int = 1 - AT1_preaff_scale/2;
AT1_effect_on_preaff = AT1_preaff_int + AT1_preaff_scale/(1+exp(-(AT1_bound_AngII - nominal_equilibrium_AT1_bound_AngII)/AT1_preaff_slope));

AT1_aff_int = 1 - AT1_aff_scale/2;
AT1_effect_on_aff = AT1_aff_int + AT1_aff_scale/(1+exp(-(AT1_bound_AngII - nominal_equilibrium_AT1_bound_AngII)/AT1_aff_slope));

AT1_eff_int = 1 - AT1_eff_scale/2;
AT1_effect_on_eff = AT1_eff_int + AT1_eff_scale/(1+exp(-(AT1_bound_AngII - nominal_equilibrium_AT1_bound_AngII)/AT1_eff_slope));


# ANP_aff_int = 1 + ANP_aff_scale/2;
# ANP_effect_on_aff = ANP_aff_int - ANP_aff_scale/(1+exp(-(ANP - nom_ANP)/ANP_aff_slope));

### RSNA constricts the preafferent vasculature
rsna_preaff_int = 1 - rsna_preaff_scale/2;
rsna_effect_on_preaff = rsna_preaff_int + rsna_preaff_scale/(1+exp(-(renal_sympathetic_nerve_activity - nom_rsna)/rsna_preaff_slope));

### ANP may dilate the preafferent, afferent, and/or efferent arterioles
# ANP_preaff_int = 1 + ANP_preaff_scale/2;
# ANP_effect_on_preaff = ANP_preaff_int - ANP_preaff_scale/(1+exp(-(normalized_atrial_NP_concentration - nom_ANP)/ANP_preaff_slope));

# ANP_aff_int = 1 + ANP_aff_scale/2;
# ANP_effect_on_aff = ANP_aff_int - ANP_aff_scale/(1+exp(-(normalized_atrial_NP_concentration - nom_ANP)/ANP_aff_slope));

# ANP_eff_int = 1 + ANP_eff_scale/2;
# ANP_effect_on_eff= ANP_eff_int - ANP_eff_scale/(1+exp(-(normalized_atrial_NP_concentration - nom_ANP)/ANP_eff_slope));


######Preafferent Resistance


if(renal_blood_flow_L_min_delayed < 1 & interstitial_fluid_volume < 15){
IF_Venous_RIHP_Effect_int = 1- IF_Venous_RIHP_Effect_scale/2;
IF_Venous_RIHP_Effect =(IF_Venous_RIHP_Effect_int+IF_Venous_RIHP_Effect_scale/(1+exp(-(renal_blood_flow_L_min_delayed-nom_renal_blood_flow_L_min)/IF_Venous_RIHP_Effect_slope)));
}else{
IF_Venous_RIHP_Effect=1;
}


preaff_arteriole_signal_multiplier = AT1_effect_on_preaff*rsna_effect_on_preaff*preafferent_pressure_autoreg_signal*CCB_effect_on_preafferent_resistance;
preaff_arteriole_adjusted_signal_multiplier = (1/(1+exp(preaff_signal_nonlin_scale*(1-preaff_arteriole_signal_multiplier)))+0.5);
preafferent_arteriole_resistance = IF_Venous_RIHP_Effect*nom_preafferent_arteriole_resistance*preaff_arteriole_adjusted_signal_multiplier;

###### Afferent Arteriole Resistance
#The afferent arteriole responses the tubuloglomerular feedback (calculated later), as well as to AT1-bound AngII and ANP.
#It may respond myogenically as well. Some studies suggest the upstream portion responds myogenically while the distal portion responds to TGF. Thus, one could consider the
#myogenically responsive portion as part of the preafferent resistance.
#The dilation/constriction of the arterioles is limited, and thus the total combined effect of all regulators must saturate

# afferent_arteriole_signal_multiplier = tubulo_glomerular_feedback_effect * AT1_effect_on_aff * ANP_effect_on_aff*glomerular_pressure_autoreg_signal*CCB_effect_on_afferent_resistance;
afferent_arteriole_signal_multiplier = tubulo_glomerular_feedback_effect * AT1_effect_on_aff *glomerular_pressure_autoreg_signal*CCB_effect_on_afferent_resistance;
afferent_arteriole_adjusted_signal_multiplier = (1/(1+exp(afferent_signal_nonlin_scale*(1-afferent_arteriole_signal_multiplier)))+0.5);
afferent_arteriole_resistance = IF_Venous_RIHP_Effect*nom_afferent_arteriole_resistance*afferent_arteriole_adjusted_signal_multiplier*(1-ANP_effect_on_Arterial_Resistance);

###### Efferent Arteriole Resistance
#The efferent arteriole responses to AT1-bound AngII and ANP.
#The dilation/constriction of the arterioles is limited, and thus the total combined effect of all regulators must saturate
#efferent_arteriole_signal_multiplier = AT1_effect_on_eff * ANP_effect_on_eff *CCB_effect_on_efferent_resistance;

efferent_arteriole_signal_multiplier = AT1_effect_on_eff * CCB_effect_on_efferent_resistance;
efferent_arteriole_adjusted_signal_multiplier = 1/(1+exp(efferent_signal_nonlin_scale*(1-efferent_arteriole_signal_multiplier)))+0.5;
efferent_arteriole_resistance =  IF_Venous_RIHP_Effect*nom_efferent_arteriole_resistance*efferent_arteriole_adjusted_signal_multiplier; #

######Peritubular Resistance
#Autoregulation of peritubular resistance allows RBF to be autoregulated separately from GFR
#This is exploratory for now. By default, this effect is turned off by setting RBF_autoreg_scale to zero
RBF_autoreg_int = 1-RBF_autoreg_scale/2;
peritubular_autoreg_signal = RBF_autoreg_int + RBF_autoreg_scale/(1+exp((nom_renal_blood_flow_L_min - renal_blood_flow_L_min_delayed)/RBF_autoreg_steepness));
autoregulated_peritubular_resistance = peritubular_autoreg_signal*nom_peritubular_resistance;

###### Renal Vascular Resistance
renal_vascular_resistance = (preafferent_arteriole_resistance + (afferent_arteriole_resistance + efferent_arteriole_resistance) / number_of_functional_glomeruli + autoregulated_peritubular_resistance);

###Renal blood flow
renal_blood_flow_L_min = ((mean_arterial_pressure_MAP - P_venous) / renal_vascular_resistance);
renal_blood_flow_ml_hr = renal_blood_flow_L_min * 1000 * 60;


###Renal Vasculature Pressures
preafferent_pressure = mean_arterial_pressure_MAP - renal_blood_flow_L_min*preafferent_arteriole_resistance;
glomerular_pressure = (mean_arterial_pressure_MAP  - renal_blood_flow_L_min * (preafferent_arteriole_resistance + afferent_arteriole_resistance / number_of_functional_glomeruli));
postglomerular_pressure = (mean_arterial_pressure_MAP  - renal_blood_flow_L_min * (preafferent_arteriole_resistance + (afferent_arteriole_resistance+efferent_arteriole_resistance) / number_of_functional_glomeruli));


#Autoregulatory signals for preafferent and afferent resistances
preaff_autoreg_int = 1 - preaff_autoreg_scale/2;
preafferent_pressure_autoreg_function = preaff_autoreg_int+preaff_autoreg_scale/(1+exp((nom_preafferent_pressure - preafferent_pressure)/myogenic_steepness));
gp_autoreg_int = 1 - gp_autoreg_scale/2;
glomerular_pressure_autoreg_function = gp_autoreg_int+gp_autoreg_scale/(1+exp((nom_glomerular_pressure - glomerular_pressure)/myogenic_steepness));

######################## Glomerular Filtration ######################

# Assume glomerulosclerosis causes a decrease in Kf over time, and also a loss of the renal vasculature (afferent and efferent arterioles)
# as glomeruli become completely sclerotic

#Glomerular hypertrophy resulting in increased surface area and thus increased Kf is assumed to occur
#in response to elevated glomerular pressure. A 2 mmHg buffer is built in (i.e. glomerular pressure must be at least 2 mmHg above normal for hypertrophy to begin
#The increase in Kf saturates and cannot exceed the fractional increase set by maximal_glom_surface_area_increase
GP_effect_increasing_Kf = (maximal_glom_surface_area_increase - disease_effects_increasing_Kf) * max(glomerular_pressure/(nom_glomerular_pressure+2) - 1,0) / (T_glomerular_pressure_increases_Kf/C_renal_CV_timescale);
glomerular_hydrostatic_conductance_Kf = nom_Kf*(1+disease_effects_increasing_Kf)*(1-disease_effects_decreasing_Kf);

###Glomerular Fitlration Rate
#Calculation of P_bowmans are described later

net_filtration_pressure = glomerular_pressure - oncotic_pressure_difference - P_bowmans;

if (net_filtration_pressure <= 0) { #Prevent crashing if NFP becomes negative
SNGFR_nL_min = 0.001;
} else {
SNGFR_nL_min = glomerular_hydrostatic_conductance_Kf * (glomerular_pressure - oncotic_pressure_difference - P_bowmans);
}

GFR =  (SNGFR_nL_min / 1000 / 1000000 * number_of_functional_tubules);
GFR_ml_min = GFR * 1000;

filtration_fraction = GFR/renal_blood_flow_L_min;

serum_creatinine_concentration = serum_creatinine/blood_volume_L;
creatinine_clearance_rate = GFR_ml_min * dl_ml * serum_creatinine_concentration; #Units: mg/min

####################### Protein filtering ####################
GPdiff = max(0, glomerular_pressure - (nom_GP_seiving_damage));
GP_effect_on_Seiving = Emax_seiving * GPdiff ^ Gamma_seiving / (GPdiff ^ Gamma_seiving + Km_seiving ^ Gamma_seiving);

#Dean and Lazzara 2006 - Seiving coefficient decreases as GFR increases

nom_glomerular_albumin_sieving_coefficient = seiving_inf/(1-(1-seiving_inf)*exp(-c_albumin*SNGFR_nL_min));
glomerular_albumin_sieving_coefficient = nom_glomerular_albumin_sieving_coefficient*(1 + GP_effect_on_Seiving);

SN_albumin_filtration_rate = plasma_albumin_concentration*SNGFR_nL_min*1e-6*glomerular_albumin_sieving_coefficient; #mg/min

SN_albumin_excretion_rate = max(0, SN_albumin_filtration_rate - SN_albumin_reabsorptive_capacity)+nom_albumin_excretion_rate;
albumin_excretion_rate = SN_albumin_excretion_rate*number_of_functional_tubules;

#albumin_filtation_rate = plasma_albumin_concentration*SNGFR_nL_min*glomerular_albumin_sieving_coefficient;
#albumin_excretion_rate = max(0, albumin_filtation_rate - max_PT_albumin_reabsorption_rate);



###Oncotic pressure difference
#Landis Pappenheimer equation used to calculate oncotic pressure at entrance and exit to glomerulus
#Oncotic pressure is approximated as varying linearly along the glomerulus. Oncotic pressure in the Bowman's space is zero
#Thus the average pressure difference is the average of the entrance and exit oncotic pressure
#We do not consider filtration equilibrium
Oncotic_pressure_in = 1.629*plasma_protein_concentration+0.2935*(plasma_protein_concentration^2);
SNRBF_nl_min = 1e6*1000*renal_blood_flow_L_min/number_of_functional_glomeruli;
plasma_protein_concentration_out = (SNRBF_nl_min*plasma_protein_concentration-SN_albumin_filtration_rate)/(SNRBF_nl_min-SNGFR_nL_min);
Oncotic_pressure_out = 1.629*plasma_protein_concentration_out+0.2935*(plasma_protein_concentration_out^2);
oncotic_pressure_avg = (Oncotic_pressure_in+Oncotic_pressure_out)/2;



##################### Plasma sodium concentration and vasopressin secretion
###Plasma sodium concentration
Na_concentration = sodium_amount / blood_volume_L;
IF_Na_concentration = IF_sodium_amount/interstitial_fluid_volume;

sodium_storate_rate = Q_Na_store*((max_stored_sodium - stored_sodium)/max_stored_sodium)*(IF_Na_concentration - ref_Na_concentration);

###Control of vasopressin secretion
#A proportional-integral controller is used to ensure there is no steady state error in sodium concentration
#Relative gains of the P and I controller must be chosen carefully.
#In order to permit a steady-state error, the integral controller can be removed. But care should be given then in choosing the proportional gain
Na_water_controller = Na_controller_gain*(Kp_VP*(Na_concentration - ref_Na_concentration)+Ki_VP*Na_concentration_error);

###Vasopressin
#Vasopressin is critical in the model, because it allows water excretion to be decoupled from sodium excretion in the collecting duct
normalized_vasopressin_concentration = 1 + Na_water_controller;
vasopressin_concentration = nominal_vasopressin_conc * normalized_vasopressin_concentration;

#Effect of vasopressin on water intake
water_intake_vasopressin_int = 1-water_intake_vasopressin_scale/2;
water_intake = water_intake_species_scale*(nom_water_intake/60/24)*(water_intake_vasopressin_int + water_intake_vasopressin_scale/(1+exp((normalized_vasopressin_concentration_delayed-1)/water_intake_vasopressin_slope)));
daily_water_intake = (water_intake * 24 * 60);


##################### Tubular Flow and Reabsorption######################

#Length of tubular segments
L_pt_s1 = L_pt_s1_nom*(1+tubular_length_increase);
L_pt_s2 = L_pt_s2_nom*(1+tubular_length_increase);
L_pt_s3 = L_pt_s3_nom*(1+tubular_length_increase);
Dc_pt = Dc_pt_nom*(1+tubular_diameter_increase);
L_pt = L_pt_s1+L_pt_s2 + L_pt_s3;

SN_filtered_Na_load = (SNGFR_nL_min / 1000 / 1000000)*Na_concentration;
filtered_Na_load = SN_filtered_Na_load*number_of_functional_tubules;

######Regulatory effects on reabsorption
###Pressure natriuresis effects


pressure_natriuresis_f_int = 1- pressure_natriuresis_f_scale/2;
pressure_natriuresis_p_int = 1- pressure_natriuresis_p_scale/2;

pressure_natriuresis_signal_f = ((pressure_natriuresis_f_int+pressure_natriuresis_f_scale/(1+exp(-(renal_blood_flow_L_min_delayed-nom_renal_blood_flow_L_min)/pressure_natriuresis_f_slope))));  #+0.85 +0.3

pressure_natriuresis_signal_p =(pressure_natriuresis_p_int+pressure_natriuresis_p_scale/(1+exp(-(RIHP_delayed-RIHP0)/pressure_natriuresis_p_slope)));

pressure_natriuresis_IF_int = 1- pressure_natriuresis_IF_scale/2;
pessure_natriuresis_signal_IF = ((pressure_natriuresis_IF_int+pressure_natriuresis_IF_scale/(1+exp(-(Net_oncotic_pressure_diff-Nom_Net_oncotic_pressure_diff)/pressure_natriuresis_IF_slope))));


pressure_natriuresis_signal = pressure_natriuresis_signal_f*pressure_natriuresis_signal_p*pessure_natriuresis_signal_IF;

pressure_natriuresis_PT_int = 1 - pressure_natriuresis_PT_scale/2;
pressure_natriuresis_PT_effect = max(0.001,pressure_natriuresis_PT_int + pressure_natriuresis_PT_scale / (1 + exp(pressure_natriuresis_signal-1)));

pressure_natriuresis_LoH_int = 1 - pressure_natriuresis_LoH_scale/2;
pressure_natriuresis_LoH_effect = max(0.001,pressure_natriuresis_LoH_int + pressure_natriuresis_LoH_scale / (1 + exp((postglomerular_pressure_delayed - RIHP0) / pressure_natriuresis_LoH_slope)));

pressure_natriuresis_DCT_magnitude = max(0,pressure_natriuresis_DCT_scale );
pressure_natriuresis_DCT_int = 1 - pressure_natriuresis_DCT_magnitude/2;
pressure_natriuresis_DCT_effect = max(0.001,pressure_natriuresis_DCT_int + pressure_natriuresis_DCT_magnitude/ (1 + exp((postglomerular_pressure_delayed - RIHP0) / pressure_natriuresis_DCT_slope)));

pressure_natriuresis_CD_magnitude = max(0,pressure_natriuresis_CD_scale *(1+disease_effects_decreasing_CD_PN));
pressure_natriuresis_CD_int = 1 - pressure_natriuresis_CD_magnitude/2;
pressure_natriuresis_CD_effect = max(0.001,pressure_natriuresis_CD_int + pressure_natriuresis_CD_magnitude/ (1 + exp(pressure_natriuresis_signal-1)));

#smax(0.001,pressure_natriuresis_signal);

###AT1-bound AngII effect on PT reabsorption
AT1_PT_int = 1 - AT1_PT_scale/2;
AT1_effect_on_PT = AT1_PT_int + AT1_PT_scale/(1+exp(-(AT1_bound_AngII - nominal_equilibrium_AT1_bound_AngII)/AT1_PT_slope));
### Aldosterone effect on DCT and CD reabsorption



### RSNA effect on PT and CD sodium reabsorption
#RSNA effect on CD is turned off by default
rsna_PT_int = 1 - rsna_PT_scale/2;

#######NEED To either change rsna_delayed consistently throughout or revert back
rsna_effect_on_PT = 1;
rsna_CD_int = 1 - rsna_CD_scale/2;
rsna_effect_on_CD= rsna_CD_int + rsna_CD_scale/(1+exp((1 - renal_sympathetic_nerve_activity)/rsna_CD_slope));


###Aldosterone effect on distal and collecting duct sodium reabsorption
aldosterone_concentration = normalized_aldosterone_level* nominal_aldosterone_concentration;
Aldo_MR_normalised_effect = normalized_aldosterone_level*(1 - pct_target_inhibition_MRA);

aldo_DCT_int = 1 - aldo_DCT_scale/2;
aldo_effect_on_DCT = aldo_DCT_int + aldo_DCT_scale/(1+exp((1 - Aldo_MR_normalised_effect)/aldo_DCT_slope));

aldo_CD_int = 1 - aldo_CD_scale/2;
aldo_effect_on_CD= aldo_CD_int + aldo_CD_scale/(1+exp((1 - Aldo_MR_normalised_effect)/aldo_CD_slope));

###ANP effect on collecting duct sodium reabsorption
# anp_CD_int = 1 - anp_CD_scale/2;
# anp_effect_on_CD= anp_CD_int + anp_CD_scale/(1+exp((1 - normalized_atrial_NP_concentration)/anp_CD_slope));

#Assume insulin has effect on NHE3. Use RUGE as surrogate for insulin effect. When RUGE goes up, insulin effect goes down.
NHE3inhib = Anhe3*RUGE_delayed;
pt_multiplier = AT1_effect_on_PT * rsna_effect_on_PT *pressure_natriuresis_PT_effect*(1-NHE3inhib);

e_pt_sodreab = min(1,nominal_pt_na_reabsorption_nonSGLT * pt_multiplier);# AT1_effect_on_PT * rsna_effect_on_PT *pressure_natriuresis_PT_effect;


e_dct_sodreab = min(1,nominal_dt_na_reabsorption * aldo_effect_on_DCT*pressure_natriuresis_DCT_effect ); #HCTZ_effect_on_DT_Na_reabs

# cd_multiplier = anp_effect_on_CD *aldo_effect_on_CD*rsna_effect_on_CD*pressure_natriuresis_CD_effect;
cd_multiplier = aldo_effect_on_CD*rsna_effect_on_CD*pressure_natriuresis_CD_effect;
#Ensure that CD reabsorption smoothly approaches upper limit set by max_cd_reabs_rate
cd_scale = max_cd_reabs_rate/nominal_cd_na_reabsorption-1;
#if (cd_multiplier > 1) {
#	cd_multiplier = 1+cd_scale - cd_scale*exp((1-cd_multiplier)/.028);
#}
e_cd_sodreab = min(0.9999,nominal_cd_na_reabsorption*cd_multiplier*lambdaANP);

##################################### Proximal Tubule #########################################
###Glucose Filtration and reabsorption in PT
#Assume glucose reabsorption depends only on availability of SGLT1/2
#Assume constant amount of reabsorption per unit length through SGLT2 in convoluted PT
#Assume constant amount of reabsorption per unit length through SGLT1 in straight/recta PT

#Chosen so that UGE becomes non-zero for plasma_glucose concentration ~8.5 mmol/l
glucose_reabs_per_unit_length_s1 = nom_glucose_reabs_per_unit_length_s1*diabetic_adaptation*SGLT2_inhibition_delayed*(1+RTg_compensation);
glucose_reabs_per_unit_length_s2 = nom_glucose_reabs_per_unit_length_s2*diabetic_adaptation*SGLT2_inhibition_delayed*(1+RTg_compensation);
glucose_reabs_per_unit_length_s3 = nom_glucose_reabs_per_unit_length_s3*diabetic_adaptation*(1+RTg_compensation)*SGLT1_inhibition;

SN_filtered_glucose_load = glucose_concentration*SNGFR_nL_min / 1000 / 1000000;  #mmol/min
glucose_pt_out_s1 = max(0,SN_filtered_glucose_load-glucose_reabs_per_unit_length_s1*L_pt_s1); #mmol/min
glucose_pt_out_s2 = max(0,glucose_pt_out_s1-glucose_reabs_per_unit_length_s2*L_pt_s2); #mmol/min
glucose_pt_out_s3 = max(0,glucose_pt_out_s2-glucose_reabs_per_unit_length_s3*L_pt_s3); #mmol/min
RUGE = glucose_pt_out_s3*number_of_functional_tubules*180; #RUGE in mg/min


excess_glucose_increasing_RTg = (maximal_RTg_increase - RTg_compensation) * max(RUGE,0) / (T_glucose_RTg/C_renal_CV_timescale);

osmotic_natriuresis_effect_pt = 1-min(1,RUGE *glucose_natriuresis_effect_pt);
osmotic_natriuresis_effect_cd = 1-min(1,RUGE *glucose_natriuresis_effect_cd);
osmotic_diuresis_effect_pt = 1-min(1,RUGE *glucose_diuresis_effect_pt);
osmotic_diuresis_effect_cd = 1-min(1,RUGE *glucose_diuresis_effect_cd);


###PT Sodium filtration and reabsorption
# Sodium reabsorbed 1:1 with glucose in S1 and S2
# Sodium reabsorbed 2:1 with glucose in S3
# Assume for non-SGLT reabsorption, sodium reabsorbed at a constant RATE along the tubule
# (represents glomerulotubular balance)
SN_filtered_Na_load = (SNGFR_nL_min / 1000 / 1000000)*Na_concentration;  	#mmol/min
SGTL2_Na_reabs_mmol_s1 = SN_filtered_glucose_load-glucose_pt_out_s1;   		#mmol/min
SGTL2_Na_reabs_mmol_s2 = glucose_pt_out_s1-glucose_pt_out_s2;			#mmol/min
SGTL1_Na_reabs_mmol = 2*(glucose_pt_out_s2-glucose_pt_out_s3);			#mmol/min
total_SGLT_Na_reabs = SGTL2_Na_reabs_mmol_s1+SGTL2_Na_reabs_mmol_s2+SGTL1_Na_reabs_mmol; 	#mmol/min



Na_reabs_per_unit_length = -log(1-e_pt_sodreab)/(L_pt_s1+L_pt_s2+L_pt_s3); #non-SGLT2 reabs	#mmol/min
Na_pt_s1_reabs = min(max_s1_Na_reabs, SN_filtered_Na_load*(1-exp(-Na_reabs_per_unit_length*L_pt_s1)));
Na_pt_out_s1 = SN_filtered_Na_load - Na_pt_s1_reabs - SGTL2_Na_reabs_mmol_s1 ;

Na_pt_s2_reabs = min(max_s2_Na_reabs, Na_pt_out_s1*(1-exp(-Na_reabs_per_unit_length*L_pt_s2)));
Na_pt_out_s2 = Na_pt_out_s1 - Na_pt_s2_reabs - SGTL2_Na_reabs_mmol_s2;

Na_pt_s3_reabs = min(max_s3_Na_reabs, Na_pt_out_s2*(1-exp(-Na_reabs_per_unit_length*L_pt_s3)));
Na_pt_out_s3 = Na_pt_out_s2 - Na_pt_s3_reabs - SGTL1_Na_reabs_mmol;

PT_Na_reabs_fraction = 1-Na_pt_out_s3/SN_filtered_Na_load;

###PT Urea filtration and reabsorption
SN_filtered_urea_load = (SNGFR_nL_min / 1000 / 1000000)*plasma_urea;
urea_out_s1 = SN_filtered_urea_load - urea_permeability_PT*(SN_filtered_urea_load/(0.5*((SNGFR_nL_min / 1000 / 1000000)+water_out_s1_delayed))-plasma_urea)*water_out_s1_delayed; #For now, assuming only reabsorbed at the end
urea_out_s2 = urea_out_s1 - urea_permeability_PT*(urea_out_s1/(0.5*(water_out_s1_delayed+water_out_s2_delayed))-plasma_urea)*water_out_s2_delayed; #For now, assuming only reabsorbed at the end
urea_out_s3 = urea_out_s2 - urea_permeability_PT*(urea_out_s2/(0.5*(water_out_s2_delayed+water_out_s3_delayed))-plasma_urea)*water_out_s3_delayed; #For now, assuming only reabsorbed at the end
urea_reabsorption_fraction = 1-urea_out_s3/SN_filtered_urea_load;


###PT Water Reabsorption
osmoles_out_s1 = 2*Na_pt_out_s1 + glucose_pt_out_s1 + urea_out_s1;
water_out_s1 = (((SNGFR_nL_min / 1000 / 1000000)/(2*SN_filtered_Na_load+SN_filtered_glucose_load+ SN_filtered_urea_load)))*osmoles_out_s1;

osmoles_out_s2 = 2*Na_pt_out_s2 + glucose_pt_out_s2 + urea_out_s2;
water_out_s2 = (water_out_s1/osmoles_out_s1)*osmoles_out_s2;

osmoles_out_s3 = 2*Na_pt_out_s3 + glucose_pt_out_s3 + urea_out_s3;
water_out_s3 = (water_out_s2/osmoles_out_s2)*osmoles_out_s3;

PT_water_reabs_fraction = 1-water_out_s3/(SNGFR_nL_min / 1000 / 1000000);

###Concentrations out of PT
Na_concentration_out_s1 = Na_pt_out_s1/water_out_s1;
Na_concentration_out_s2 = Na_pt_out_s2/water_out_s2;
Na_concentration_out_s3 = Na_pt_out_s3/water_out_s3;

glucose_concentration_out_s1 = glucose_pt_out_s1/water_out_s1;
glucose_concentration_out_s2 = glucose_pt_out_s2/water_out_s2;
glucose_concentration_out_s3 = glucose_pt_out_s3/water_out_s3;

urea_concentration_out_s1 = urea_out_s1/water_out_s1;
urea_concentration_out_s2 = urea_out_s2/water_out_s2;
urea_concentration_out_s3 = urea_out_s3/water_out_s3;

osmolality_out_s1 = osmoles_out_s1/water_out_s1;
osmolality_out_s2 = osmoles_out_s2/water_out_s2;
osmolality_out_s3 = osmoles_out_s3/water_out_s3;

PT_Na_outflow = Na_pt_out_s3*number_of_functional_tubules;

#Tubular sodium reabsorption per unit SA as the driver of tubular hypertrophy
PT_Na_reab_perUnitSA = SN_filtered_Na_load*e_pt_sodreab/(3.14*Dc_pt*(L_pt_s1+L_pt_s2+L_pt_s3));

normalized_PT_reabsorption_density = PT_Na_reab_perUnitSA/PT_Na_reab_perUnitSA_0;
PT_Na_reabs_effect_increasing_tubular_length = 0;#(maximal_tubule_length_increase - tubular_length_increase) * max(normalized_PT_reabsorption_density - 1,0) / (T_PT_Na_reabs_PT_length/C_renal_CV_timescale);
PT_Na_reabs_effect_increasing_tubular_diameter = 0;#(maximal_tubule_diameter_increase - tubular_diameter_increase) * max(normalized_PT_reabsorption_density - 1,0) / (T_PT_Na_reabs_PT_diameter/C_renal_CV_timescale);


##################################### Loop of Henle #########################################

#####Descending Loop of Henle

water_in_DescLoH = water_out_s3; # L/min
Na_in_DescLoH = Na_pt_out_s3;
urea_in_DescLoH = urea_out_s3;
glucose_in_DescLoH = glucose_pt_out_s3;
osmoles_in_DescLoH = osmoles_out_s3;

Na_concentration_in_DescLoH = Na_concentration_out_s3;
Urea_concentration_in_DescLoH = urea_concentration_out_s3;
glucose_concentration_in_DescLoH = glucose_concentration_out_s3;
osmolality_in_DescLoH = osmoles_out_s3/water_out_s3;

#No solute reabsorption in descending limb
Na_out_DescLoH = Na_in_DescLoH;
urea_out_DescLoH = urea_in_DescLoH;
glucose_out_DescLoH = glucose_in_DescLoH;
osmoles_out_DescLoH = osmoles_in_DescLoH;


#For LoH, baseline osmoles reabsorbed per unit length is calculated from nominal fractional sodium reabsorption (see baseline parameters file)
#The rate of reabsorption per unit length may be flow-dependent, and may be modulated by tubular pressure-natriuresis
# If LoH_flow_dependence = 0, then no flow dependence.

deltaLoH_NaFlow = min(max_deltaLoH_reabs,LoH_flow_dependence*(Na_out_DescLoH-nom_Na_in_AscLoH));
AscLoH_Reab_Rate =(2*nominal_loh_na_reabsorption*(nom_Na_in_AscLoH+deltaLoH_NaFlow)*loop_diuretic_effect)/L_lh_des; #osmoles reabsorbed per unit length per minute. factor of 2 because osmoles = 2

effective_AscLoH_Reab_Rate =AscLoH_Reab_Rate*pressure_natriuresis_LoH_effect; #osmoles reabsorbed per unit length per minute. factor of 2 because osmoles = 2*Na


#Min function necesssary to ensure that the LoH does not reabsorb more Na than is delivered to it
osmolality_out_DescLoH = osmolality_in_DescLoH*exp(min(effective_AscLoH_Reab_Rate*L_lh_des,2*Na_in_DescLoH)/(water_in_DescLoH*osmolality_in_DescLoH));
water_out_DescLoH = water_in_DescLoH*osmolality_in_DescLoH/osmolality_out_DescLoH;

Na_concentration_out_DescLoH = Na_out_DescLoH/water_out_DescLoH;
glucose_concentration_out_DescLoH = glucose_out_DescLoH/water_out_DescLoH;
urea_concentration_out_DescLoH = urea_out_DescLoH/water_out_DescLoH;

#####Ascending Loop of Henle

Na_in_AscLoH = Na_out_DescLoH;
urea_in_AscLoH_before_secretion = urea_out_DescLoH;
glucose_in_AscLoH = glucose_out_DescLoH;
osmoles_in_AscLoH_before_secretion = osmoles_out_DescLoH;
water_in_AscLoH = water_out_DescLoH;
Na_concentration_in_AscLoH = Na_concentration_out_DescLoH;
#Urea Secretion --> Assume all urea reabsorbed and secreted only at tip of loop
urea_in_AscLoH = urea_in_AscLoH_before_secretion + reabsorbed_urea_cd_delayed;
urea_concentration_in_AscLoH = urea_in_AscLoH/water_out_DescLoH;
osmoles_in_AscLoH = osmoles_in_AscLoH_before_secretion  + reabsorbed_urea_cd_delayed;

osmolality_in_AscLoH = osmoles_in_AscLoH/water_in_AscLoH;

#Osmolality descreased due to sodium reabsorption along ascending loop
#min function necessary so that LoH doesn't reabsorb more sodium than is delivered to it
osmolality_out_AscLoH = osmolality_in_AscLoH - min(L_lh_des*effective_AscLoH_Reab_Rate, 2*Na_in_DescLoH)*(exp(min(L_lh_des*effective_AscLoH_Reab_Rate, 2*Na_in_DescLoH)/(water_in_DescLoH*osmolality_in_DescLoH))/water_in_DescLoH);
osmoles_reabsorbed_AscLoH = (osmolality_in_AscLoH - osmolality_out_AscLoH)*water_in_AscLoH;
Na_reabsorbed_AscLoH = osmoles_reabsorbed_AscLoH/2;
Na_out_AscLoH = max(0,Na_in_AscLoH - Na_reabsorbed_AscLoH);



#Water, glucose, and urea are not reabsorbed along the ascending limb
urea_out_AscLoH = urea_in_AscLoH; #urea secretion accounted for above
glucose_out_AscLoH = glucose_in_AscLoH;
water_out_AscLoH = water_in_AscLoH;

osmoles_out_AscLoH = osmolality_out_AscLoH*water_out_AscLoH;


Na_concentration_out_AscLoH = Na_out_AscLoH/water_out_AscLoH;
glucose_concentration_out_AscLoH = glucose_out_AscLoH/water_out_AscLoH;
urea_concentration_out_AscLoH = urea_out_AscLoH/water_out_AscLoH;

LoH_reabs_fraction = 1-Na_out_AscLoH/Na_in_AscLoH;

SN_macula_densa_Na_flow = Na_out_AscLoH;
MD_Na_concentration = Na_concentration_out_AscLoH;
TGF0_tubulo_glomerular_feedback = 1 - S_tubulo_glomerular_feedback/2;
tubulo_glomerular_feedback_signal = (TGF0_tubulo_glomerular_feedback + S_tubulo_glomerular_feedback / (1 + exp((MD_Na_concentration_setpoint - MD_Na_concentration)/ F_md_scale_tubulo_glomerular_feedback)));


#############################Distal Convoluted Tubule #######################


water_in_DCT = water_out_AscLoH;
Na_in_DCT = Na_out_AscLoH;
urea_in_DCT = urea_out_AscLoH;
glucose_in_DCT = glucose_out_AscLoH;
osmoles_in_DCT = osmoles_out_AscLoH;

Na_concentration_in_DCT = Na_concentration_out_AscLoH;
urea_concentration_in_DCT = urea_concentration_out_AscLoH;
glucose_concentration_in_DCT = glucose_concentration_out_AscLoH;
osmolality_in_DCT = osmolality_out_AscLoH;

#Assume only sodium reabsorbed along DCT, no water, glucose, or urea reabsorption
urea_out_DCT = urea_in_DCT;
glucose_out_DCT = glucose_in_DCT;
water_out_DCT = water_in_DCT;

urea_concentration_out_DCT = urea_out_DCT/water_out_DCT;
glucose_concentration_out_DCT = glucose_out_DCT/water_out_DCT;

#Assume sodium reabsorption at a constant fraction of delivery
R_dct = -log(1-e_dct_sodreab)/L_dct;

Na_out_DCT = Na_in_DCT*exp(-R_dct*L_dct);
Na_concentration_out_DCT = Na_out_DCT/water_out_DCT;

osmolality_out_DCT = 2*Na_concentration_out_DCT + glucose_concentration_out_DescLoH + urea_concentration_in_AscLoH;
osmoles_out_DCT = osmolality_out_DCT*water_out_DCT;

DCT_Na_reabs_fraction = 1-Na_out_DCT/Na_in_DCT;

Na_reabsorbed_DCT = Na_in_DCT - Na_out_DCT;




################################################Collecting Duct###############################

water_in_CD = water_out_DCT;
Na_in_CD = Na_out_DCT;
urea_in_CD = urea_out_DCT;
glucose_in_CD = glucose_out_DCT;

osmoles_in_CD = osmoles_out_DCT;

#Use this to turn off osmotic diuresis effect
#osmoles_in_CD = osmoles_out_DCT - glucose_in_CD;

osmolality_in_CD = osmoles_in_CD/water_in_CD;

Na_concentration_in_CD = Na_concentration_out_DCT;
urea_concentration_in_CD = urea_concentration_out_DCT;
glucose_concentration_in_CD = glucose_concentration_out_DCT;

osmotic_diuresis_effect_cd = 1-min(1,RUGE *glucose_diuresis_effect_cd);

####Assume sodium reabsorbed, then water follows
#### Then urea reabsorbed at end
#### Then additional water reabsorbed following urea reabsorption

#Assume sodium reabsorbed at fractional rate eta
e_cd_sodreab_adj = e_cd_sodreab*osmotic_natriuresis_effect_cd;
R_cd = -log(1-e_cd_sodreab_adj)/L_cd;
Na_reabsorbed_CD = min(Na_in_CD*(1-exp(-R_cd*L_cd)),CD_Na_reabs_threshold);
Na_out_CD = Na_in_CD-Na_reabsorbed_CD;

CD_Na_reabs_fraction = 1-Na_out_CD/Na_in_CD;

Na_reabsorbed_CT=Na_in_CD-Na_out_CD;

ADH_water_permeability_old = min(0.99999,max(0,nom_ADH_water_permeability*normalized_vasopressin_concentration));
ADH_water_permeability = normalized_vasopressin_concentration/(0.15+normalized_vasopressin_concentration);


osmoles_out_CD = osmoles_in_CD-2*(Na_in_CD - Na_out_CD);

osmolality_out_CD_before_osmotic_reabsorption = osmoles_out_CD/water_in_CD;
water_reabsorbed_CD = ADH_water_permeability*osmotic_diuresis_effect_cd*water_in_CD*(1-osmolality_out_CD_before_osmotic_reabsorption/osmolality_out_DescLoH);
water_out_CD = water_in_CD-water_reabsorbed_CD;
Na_concentration_out_CD = Na_out_CD/water_out_CD;
osmolality_out_CD_after_osmotic_reabsorption = osmoles_out_CD/water_out_CD;
glucose_concentration_after_urea_reabsorption = glucose_in_CD/water_out_CD;

urine_flow_rate = water_out_CD*number_of_functional_tubules;

daily_urine_flow = (urine_flow_rate * 60 * 24);

Na_excretion_via_urine = Na_out_CD*number_of_functional_tubules;
Na_balance = Na_intake_rate - Na_excretion_via_urine;
water_balance = daily_water_intake - daily_urine_flow;


total_NA_reabsorbed = (total_SGLT_Na_reabs +Na_reabsorbed_AscLoH +Na_reabsorbed_DCT +Na_reabsorbed_CT  )*number_of_functional_tubules;

Na_concentration_average_PT = (Na_concentration_out_s1 + Na_concentration_out_s2+Na_concentration_out_s3)/3;
Na_concentration_average_DescLoH = 0.5*(Na_concentration_in_DescLoH + Na_concentration_out_DescLoH);
Na_concentration_average_AscLoH = 0.5*(Na_concentration_in_AscLoH + Na_concentration_out_AscLoH );
Na_concentration_average_DCT = 0.5*(Na_concentration_in_DCT + Na_concentration_out_DCT);
Na_concentration_average_CD = 0.5*(Na_concentration_in_CD + Na_concentration_out_CD);

Na_concentration_average_tubule = 0.2*(Na_concentration_average_PT+Na_concentration_average_DescLoH+Na_concentration_average_AscLoH+Na_concentration_average_DCT+Na_concentration_average_CD);

Oncotic_pressure_tubule = Na_concentration_average_tubule*19.3*2;

Na_amount_in_renal_interstitium = total_NA_reabsorbed;
Na_concentration_renal_interstitium = Na_amount_in_renal_interstitium/(0.03831672); #Total kidney volume 202ml. Interstitium = 0.13*KV
Oncotic_pressure_renal_interstitium = Na_concentration_renal_interstitium*19.3*2;

Net_oncotic_pressure = Oncotic_pressure_tubule-Oncotic_pressure_renal_interstitium;

FENA = Na_excretion_via_urine/filtered_Na_load;

PT_fractional_glucose_reabs = (SN_filtered_glucose_load - glucose_pt_out_s3)/SN_filtered_glucose_load;
PT_fractional_Na_reabs = (SN_filtered_Na_load - Na_pt_out_s3)/SN_filtered_Na_load;
PT_fractional_urea_reabs = ( SN_filtered_urea_load - urea_out_s3)/SN_filtered_urea_load;
PT_fractional_water_reabs = ((SNGFR_nL_min / 1000 / 1000000) - water_out_s3)/(SNGFR_nL_min / 1000 / 1000000);

LoH_fractional_Na_reabs = (Na_in_DescLoH - Na_out_AscLoH)/Na_in_DescLoH;
LoH_fractional_urea_reabs = (urea_in_DescLoH-urea_out_AscLoH)/urea_in_DescLoH;
LoH_fractional_water_reabs = (water_in_DescLoH - water_out_AscLoH)/water_in_DescLoH;

DCT_fractional_Na_reabs = (Na_in_DCT - Na_out_DCT)/Na_in_DCT;

CD_fractional_Na_reabs = (Na_in_CD - Na_out_CD)/Na_in_CD;
#CD_fractional_urea_reabs = (urea_in_CD - urea_out_CD)/urea_in_CD;
CD_OM_fractional_water_reabs = (water_in_CD - water_out_CD)/water_in_CD;


#####################Renal Interstitial Hydrostatic pressure

######RIHP can be approximated from Starling's equation for the peritubular capillaries
### Flow out of the capillary = Kf_peritubular*(Peritubular pressure - RIHP - oncotic pressure difference)
### Assume that any fluid reabsorbed reenters the capillary.
### Therefore, RIHP = Peritubular Pressure - (oncotic pressure in peritubular capillary - interstitial oncotic pressure) + tubular reabsorption/KF
#Peritubular pressure is assumed to equal postglomerular pressure

#Oncotic pressure at the entrance to the peritubular circulation equals oncotic pressure at the exit of the glomerular
Oncotic_pressure_peritubular_in = Oncotic_pressure_out;
plasma_protein_concentration_peritubular_out = (SNRBF_nl_min)*plasma_protein_concentration/(SNRBF_nl_min-urine_flow_rate*1e6*1000/number_of_functional_glomeruli);
Oncotic_pressure_peritubular_out = 1.629*plasma_protein_concentration_peritubular_out+0.2935*(plasma_protein_concentration_peritubular_out^2);
oncotic_pressure_peritubular_avg = (Oncotic_pressure_peritubular_in+Oncotic_pressure_peritubular_out)/2;

Na_concentration_peritubular_cap = (sodium_amount-filtered_Na_load)/blood_volume_L;
oncotic_pressure_peritubular_cap_Na = 0;#Na_concentration_peritubular_cap*19.3*2;
oncotic_pressure_peritubular = oncotic_pressure_peritubular_avg+oncotic_pressure_peritubular_cap_Na;

tubular_reabsorption = GFR_ml_min/1000 - urine_flow_rate;

Volume = interstitial_fluid_volume;

volume_RIHP_int = 1- volume_RIHP_scale/2;
IF_effect_RIHP = (volume_RIHP_int+volume_RIHP_scale/(1+exp(-(Volume-IF_nom)/volume_RIHP_slope)));

oncotic_int = 1-IF_Interonvotic_Effect_scale/2;
IF_effect_oncotic = (oncotic_int+IF_Interonvotic_Effect_scale/(1+exp((Volume-IF_nom)/IF_Interonvotic_Effect_slope)));

#########################################################################

Renal_plasma_amount= 2.5 * RISF_nom*0.01;  # plasma amount in renal interstitium   # Plasma protein concentration = 7
RISF_plasma_protein_concentration = Renal_plasma_amount / (RISF*10);
interstitial_oncotic_pressure = 1.629*RISF_plasma_protein_concentration+0.2935*(RISF_plasma_protein_concentration^2);
RIHP = RISF/C_RISF;
capillary_filtration = nom_peritubular_cap_Kf*(RIHP - postglomerular_pressure - (interstitial_oncotic_pressure -oncotic_pressure_peritubular));

#RIHP = IF_effect_RIHP*(postglomerular_pressure - (oncotic_pressure_peritubular - interstitial_oncotic_pressure*IF_effect_oncotic) + tubular_reabsorption/nom_peritubular_cap_Kf);


################# Tubular Pressure #####################

#####See written documentation for derivation of the equations below
#flow rates expressed in m3/min, rather than L/min


mmHg_Nperm2_conv = 133.32;
Pc_pt_s1 = Pc_pt_s1_mmHg*mmHg_Nperm2_conv;
Pc_pt_s2 = Pc_pt_s2_mmHg*mmHg_Nperm2_conv;
Pc_pt_s3 = Pc_pt_s3_mmHg*mmHg_Nperm2_conv;
Pc_lh_des = Pc_lh_des_mmHg*mmHg_Nperm2_conv;
Pc_lh_asc = Pc_lh_asc_mmHg*mmHg_Nperm2_conv;
Pc_dt = Pc_dt_mmHg*mmHg_Nperm2_conv;
Pc_cd = Pc_cd_mmHg*mmHg_Nperm2_conv;
P_interstitial = (RIHP)*mmHg_Nperm2_conv;# 4.9*mmHg_Nperm2_conv;
pi=3.14;

###CD
B1 = (4*tubular_compliance+1)*128*gamma/pi;
mean_cd_water_flow = (water_in_CD-water_out_CD)/2;
B2_cd = (Pc_cd^(4*tubular_compliance))/(Dc_cd^4);
P_in_cd = (0^(4*tubular_compliance+1)+B1*B2_cd*(mean_cd_water_flow/1e3)*L_cd)^(1/(4*tubular_compliance+1));
P_in_cd_mmHg = (P_in_cd+P_interstitial)/mmHg_Nperm2_conv;



### DCT
B2_dt = (Pc_dt^(4*tubular_compliance))/(Dc_dt^4);
P_in_dt = (P_in_cd^(4*tubular_compliance+1)+B1*B2_dt*(water_in_DCT/1e3)*L_dct)^(1/(4*tubular_compliance+1));
P_in_dt_mmHg = (P_in_dt+P_interstitial)/mmHg_Nperm2_conv;



### Asc LoH
B2_lh_asc = (Pc_lh_asc^(4*tubular_compliance))/(Dc_lh^4);
P_in_lh_asc = (P_in_dt^(4*tubular_compliance+1)+B1*B2_lh_asc*(water_in_AscLoH/1e3)*L_lh_asc)^(1/(4*tubular_compliance+1));
P_in_lh_asc_mmHg = (P_in_lh_asc+P_interstitial)/mmHg_Nperm2_conv;

### Desc LoH
A_lh_des = effective_AscLoH_Reab_Rate/(water_in_DescLoH*osmolality_in_DescLoH);
B2_lh_des = (Pc_lh_des^(4*tubular_compliance))*(water_in_DescLoH/1e3)/((Dc_lh^4)*A_lh_des);
P_in_lh_des = (P_in_lh_asc^(4*tubular_compliance+1)+B1*B2_lh_des*(1-exp(-A_lh_des*L_lh_des)))^(1/(4*tubular_compliance+1));
P_in_lh_des_mmHg = (P_in_lh_des+P_interstitial)/mmHg_Nperm2_conv;

### PT

#Treat urea as if reabsorbed linearly along whole length of PT
Rurea = (SN_filtered_urea_load - urea_out_s3)/(L_pt_s1+L_pt_s2+L_pt_s3);
urea_in_s2 = SN_filtered_urea_load - Rurea*L_pt_s1;
urea_in_s3 = SN_filtered_urea_load - Rurea*(L_pt_s1+L_pt_s2);
A_na = Na_reabs_per_unit_length;

flow_integral_s3 = 2*(Na_pt_out_s2/A_na)*(1-exp(-A_na*L_pt_s3)) - (3/2)*glucose_pt_out_s2*L_pt_s3^2 + urea_in_s3*L_pt_s3 - (1/2)*Rurea*(L_pt_s3^2);
flow_integral_s2 = 2*(Na_pt_out_s1/A_na)*(1-exp(-A_na*L_pt_s2)) - (1/2)*glucose_pt_out_s1*L_pt_s2^2 + urea_in_s2*L_pt_s2 - (1/2)*Rurea*(L_pt_s2^2);
flow_integral_s1 = 2*(SN_filtered_Na_load/A_na)*(1-exp(-A_na*L_pt_s1)) - (1/2)*SN_filtered_glucose_load*L_pt_s1^2 + SN_filtered_urea_load*L_pt_s1 - (1/2)*Rurea*(L_pt_s1^2);


#S3 segment
B2_pt_s3 = (Pc_pt_s3^(4*tubular_compliance))/(Dc_pt^4);
B3_pt_s3 = (water_out_s2/1e3)/osmoles_out_s2;
P_in_pt_s3= (P_in_lh_des^(4*tubular_compliance+1)+B1*B2_pt_s3*B3_pt_s3*flow_integral_s3)^(1/(4*tubular_compliance+1));
P_in_pt_s3_mmHg = (P_in_pt_s3+P_interstitial)/mmHg_Nperm2_conv;

B2_pt_s2 = (Pc_pt_s3^(4*tubular_compliance))/(Dc_pt^4);
B3_pt_s2 = (water_out_s1/1e3)/osmoles_out_s1;
P_in_pt_s2= (P_in_pt_s3^(4*tubular_compliance+1)+B1*B2_pt_s2*B3_pt_s2*flow_integral_s2)^(1/(4*tubular_compliance+1));
P_in_pt_s2_mmHg = (P_in_pt_s2+P_interstitial)/mmHg_Nperm2_conv;

B2_pt_s1 = (Pc_pt_s1^(4*tubular_compliance))/(Dc_pt^4);
B3_pt_s1 = (SNGFR_nL_min / 1e12)/(2*SN_filtered_Na_load+SN_filtered_glucose_load+ SN_filtered_urea_load);
P_in_pt_s1= (P_in_pt_s2^(4*tubular_compliance+1)+B1*B2_pt_s1*B3_pt_s1*flow_integral_s1)^(1/(4*tubular_compliance+1));
P_in_pt_s1_mmHg = (P_in_pt_s1+P_interstitial)/mmHg_Nperm2_conv;




####################### Aldosterone and Renin Secretion

###Aldosterone is secreted in response to AT1-bound AngII and changes in potassium or sodium concentration
#Potassium concentration is treated as a constant for now
#Empirircal relationship for Karaaslan 2005


AT1_aldo_int = 1 - AT1_aldo_slope*nominal_equilibrium_AT1_bound_AngII;
AngII_effect_on_aldo = AT1_aldo_int + AT1_aldo_slope*AT1_bound_AngII;
N_als = (K_Na_ratio_effect_on_aldo * AngII_effect_on_aldo );

###Renin is secreted in response to decreases in AT1-bound AngII and decreases in MD sodium flow

###Renin is secreted in response to decreases in AT1-bound AngII, decreases in MD sodium flow, or increases in RSNA
#RSNA effect on renin secretion
rsna_renin_intercept = 1-rsna_renin_slope;
rnsa_effect_on_renin_secretion =  rsna_renin_slope * renal_sympathetic_nerve_activity + rsna_renin_intercept;


#Macula Densa Sodium flow effect on renin secretion
#This relationship is known to be non-linear, and md_renin_tau can be calibrated based on data on changes in renin as a functoin of sodium intake
md_effect_on_renin_secretion = md_renin_A*exp(-md_renin_tau*(SN_macula_densa_Na_flow_delayed*baseline_nephrons - nom_LoH_Na_outflow));


#AT1-bound AngII feedback on renin secretion
AT1_bound_AngII_effect_on_PRA = (10 ^ (AT1_PRC_slope * log10(AT1_bound_AngII / nominal_equilibrium_AT1_bound_AngII) + AT1_PRC_yint));


#Aldo effect on renin secretion
aldo_renin_intercept = 1-aldo_renin_slope;
aldo_effect_on_renin_secretion =  aldo_renin_intercept + aldo_renin_slope*Aldo_MR_normalised_effect;

#Plasma renin activity
plasma_renin_activity = concentration_to_renin_activity_conversion_plasma* plasma_renin_concentration*(1-pct_target_inhibition_DRI);

#Renin secretion
renin_secretion_rate = (log(2)/renin_half_life)*nominal_equilibrium_PRC*AT1_bound_AngII_effect_on_PRA*md_effect_on_renin_secretion*HCTZ_effect_on_renin_secretion*aldo_effect_on_renin_secretion*BB_renin_secretion_effect;

#RAAS degradation rates
renin_degradation_rate = log(2)/renin_half_life;
AngI_degradation_rate = log(2)/AngI_half_life;
AngII_degradation_rate = log(2)/AngII_half_life;
AT1_bound_AngII_degradation_rate =  log(2)/AT1_bound_AngII_half_life;
AT2_bound_AngII_degradation_rate = log(2)/AT2_bound_AngII_half_life;

#RAAS rate constants
ACE_activity = nominal_ACE_activity*(1 - pct_target_inhibition_ACEi);
chymase_activity = nominal_chymase_activity;
AT1_receptor_binding_rate = nominal_AT1_receptor_binding_rate*(1 - pct_target_inhibition_ARB);
AT2_receptor_binding_rate = nominal_AT2_receptor_binding_rate;



Blood_volume_protein_osmotic_pressure = 1.629*plasma_protein_concentration + 0.2935*plasma_protein_concentration^2;
ISF_protein_osmotic_pressure = 1.629*ISF_protein_concentration + 0.2935*ISF_protein_concentration^2;

Blood_volume_osmotic_pressure = Blood_volume_protein_osmotic_pressure+ Na_concentration*19.3*2;
ISF_osmotic_pressire = IF_Na_concentration*19.3*2 + ISF_protein_osmotic_pressure;
Protein_sodium_filtration_pressure_grad = (mean_capillary_pressure - ISF_pressure - Blood_volume_osmotic_pressure + ISF_osmotic_pressire);

Kidney_disconnect_heart = Q_water*(Na_concentration - IF_Na_concentration);
Kidney_connect_heart = -Sodium_protein_filtration_rate_Kf*(Protein_sodium_filtration_pressure_grad)*0.001;

if (heart_renal_link == 1) {
Fluid_exchanging_function=Kidney_connect_heart;
} else{
Fluid_exchanging_function=Kidney_connect_heart;
}

peripheral_volume_change=arterial_dis_circulation_volume+capillary_circulation_volume+venules_circulation_volume;
#d/dt(HR_ratio) = heart_rate_constant*heart_rate_multiplier_adjusted;
d/dt(venous_volume) = venous_flow + C_renal_CV_timescale*(venous_volume_target - venous_volume) - tricuspid_valve_flow_rate ;
d/dt(LV_volume) = mitral_valve_flow_rate - aortic_valve_flow_rate;
d/dt(arterial_volume) = (aortic_valve_flow_rate) - (systemic_blood_flow);
#d/dt(peripheral_circulation_volume)= systemic_blood_flow - venous_flow;
d/dt(RV_volume) = (tricuspid_valve_flow_rate) - (pulmonary_valve_flow_rate);
d/dt(pulmonary_arterial_volume) = pulmonary_valve_flow_rate - pulmonary_arterial_blood_flow;
#d/dt(peripheral_circulation_volume) = (systemic_blood_flow) - venous_flow;

d/dt(arterial_dis_circulation_volume)=systemic_blood_flow - capillary_blood_flow;
d/dt(capillary_circulation_volume)=arterial_dis_blood_flow-venules_blood_flow;
d/dt(venules_circulation_volume)=capillary_blood_flow-venous_flow;

d/dt(pulmonary_venous_volume) = pulmonary_arterial_blood_flow - mitral_valve_flow_rate;
d/dt(aortic_blood_flow_delayed) = C_cycle2 * (aortic_blood_flow - aortic_blood_flow_delayed);
d/dt(pulmonary_blood_flow_delayed) = C_cycle2 * (pulmonary_blood_flow - pulmonary_blood_flow_delayed);

#################################################################################################
#################################################################################################
## Hypertrophy
d/dt(change_in_myocyte_length) = kL_hypertrophy * (LV_EDS / LV_passive_stress_along_fiber_threshhold - 1);
d/dt(change_in_myocyte_diameter) = kD_hypertrophy * (LV_active_stress_peak / LV_active_stress_threshhold - 1);
d/dt(LV_active_stress_peak) = C_cycle3 *(LV_active_stress_peak_old - LV_active_stress_peak);
d/dt(sim_time)=1;
d/dt(LV_sarcomere_length_delayed) = C_cycle* (LV_sarcomere_length - LV_sarcomere_length_delayed);
d/dt(RV_sarcomere_length_delayed) = C_cycle* (RV_sarcomere_length - RV_sarcomere_length_delayed);
d/dt(LV_EDV) = C_cycle2 * (LV_EDV_old - LV_EDV);
d/dt(LV_EDP) = C_cycle2 *(LV_EDP_old - LV_EDP);
#d/dt(LV_PSP) = C_cycle2 *(LV_PSP_old - LV_PSP);
d/dt(LV_EDS) = C_cycle2 *(LV_EDS_old - LV_EDS);

d/dt(arterial_pressure_delayed) = C_cycle2 * (arterial_pressure - arterial_pressure_delayed);
d/dt(arterial_pressure_bigger_delay) = C_cycle2 * (arterial_pressure_delayed - arterial_pressure_bigger_delay);
d/dt(systolic_pressure) = C_cycle2 * (systolic_pressure_old - systolic_pressure);
d/dt(diastolic_pressure) = C_cycle2 * (diastolic_pressure_old - diastolic_pressure);

######################computing the systolic and diastolic venous pressure##########
d/dt(venous_pressure_delayed)=C_cycle2*(venous_pressure-venous_pressure_delayed);
d/dt(venous_pressure_bigger_delay)=C_cycle2*(venous_pressure_delayed-venous_pressure_bigger_delay);
d/dt(venous_systolic_pressure)=C_cycle2*(venous_systolic_pressure_old-venous_systolic_pressure);
d/dt(venous_diastolic_pressure)=C_cycle2*(venous_diastolic_pressure_old-venous_diastolic_pressure);
d/dt(mean_venous_pressure_delayed) = 1*(mean_venous_pressure - mean_venous_pressure_delayed);

######################computing the systolic and diastolic peripheral pressure##########
d/dt(capillary_pressure_delayed)=C_cycle2*(capillary_pressure-capillary_pressure_delayed);
d/dt(capillary_pressure_bigger_delay)=C_cycle2*(capillary_pressure_delayed-capillary_pressure_bigger_delay);
d/dt(capillary_systolic_pressure)=C_cycle2*(capillary_systolic_pressure_old-capillary_systolic_pressure);
d/dt(capillary_diastolic_pressure)=C_cycle2*(capillary_diastolic_pressure_old-capillary_diastolic_pressure);
d/dt(mean_capillary_pressure_delayed) = 1*(mean_capillary_pressure - mean_capillary_pressure_delayed);


d/dt(CO) = C_co*(aortic_valve_flow_rate*60/L_m3 - CO);
d/dt(CO_delayed) = C_co_delay*(CO - CO_delayed);

#RAAS Pathway
d/dt(AngI) = plasma_renin_activity - (AngI) * (chymase_activity + ACE_activity) - (AngI) * AngI_degradation_rate;
d/dt(AngII) = AngI * (chymase_activity + ACE_activity) - AngII * AngII_degradation_rate - AngII*AT1_receptor_binding_rate - AngII* (AT2_receptor_binding_rate);
d/dt(AT1_bound_AngII) = AngII * (AT1_receptor_binding_rate) - AT1_bound_AngII_degradation_rate*AT1_bound_AngII;
d/dt(AT2_bound_AngII) = AngII * (AT2_receptor_binding_rate) - AT2_bound_AngII_degradation_rate*AT2_bound_AngII;
d/dt(plasma_renin_concentration) = renin_secretion_rate - plasma_renin_concentration * renin_degradation_rate;

#Change in Interstitial fluid volume over time is determined by the different between water intake and urine outflow
d/dt(blood_volume_L) = C_renal_CV_timescale *(water_intake- urine_flow_rate+Fluid_exchanging_function);  #
d/dt(interstitial_fluid_volume) = -C_renal_CV_timescale *Fluid_exchanging_function;

#Change in total body sodium over time is determined by the different between sodium intake and excretion
d/dt(sodium_amount) = C_renal_CV_timescale * (Na_intake_rate - Na_excretion_via_urine + Q_Na*(IF_Na_concentration - Na_concentration));
d/dt(IF_sodium_amount) = C_renal_CV_timescale *(Q_Na*(Na_concentration - IF_Na_concentration) - sodium_storate_rate);

d/dt(stored_sodium) = C_renal_CV_timescale *sodium_storate_rate;

#These equations serve only to delay the input variable by one timestep. This allows the previous value of the input variable to be used in an equation that appears
#in the code before the input variable was defined
d/dt(tubulo_glomerular_feedback_effect) =  C_renal_CV_timescale *(tubulo_glomerular_feedback_signal-tubulo_glomerular_feedback_effect);
d/dt(normalized_aldosterone_level) =  C_renal_CV_timescale *C_aldo_secretion * (N_als-normalized_aldosterone_level);
d/dt(preafferent_pressure_autoreg_signal) = C_renal_CV_timescale *100*(preafferent_pressure_autoreg_function - preafferent_pressure_autoreg_signal);
d/dt(glomerular_pressure_autoreg_signal) = 0;#C_glomerular_pressure_autoreg_signal*(glomerular_pressure_autoreg_function - glomerular_pressure_autoreg_signal);

#Error signals for PI controllers of cardiac output and sodium concentration
d/dt(CO_error) = C_renal_CV_timescale*C_co_error*(CO_delayed-CO_nom);
d/dt(Na_concentration_error) = C_renal_CV_timescale *C_Na_error*(Na_concentration - ref_Na_concentration);

#This equation allows a delay between the secretion of vasopression and its effect on water intake and tubular water reabsorption
d/dt(normalized_vasopressin_concentration_delayed)= C_renal_CV_timescale *C_vasopressin_delay*(normalized_vasopressin_concentration - normalized_vasopressin_concentration_delayed);

#TGF resetting. If C_tgf_reset = 0, no TGF resetting occurs. If it is greater than zero, the setpoint will change over time and will eventually
#come to equal the ambient MD sodium flow rate.
d/dt(F0_TGF) = C_renal_CV_timescale *C_tgf_reset*(SN_macula_densa_Na_flow*baseline_nephrons - F0_TGF);

#As above, these equations allow a variable to be used in equations that appear in the code before the variable was first defined.
d/dt(P_bowmans) = C_renal_CV_timescale *100*(P_in_pt_s1_mmHg - P_bowmans);
d/dt(oncotic_pressure_difference) = 100*(oncotic_pressure_avg - oncotic_pressure_difference);
d/dt(renal_blood_flow_L_min_delayed)=C_renal_CV_timescale*C_rbf*(renal_blood_flow_L_min - renal_blood_flow_L_min_delayed);
d/dt(SN_macula_densa_Na_flow_delayed) = C_renal_CV_timescale * C_md_flow*( SN_macula_densa_Na_flow - SN_macula_densa_Na_flow_delayed);
d/dt(rsna_delayed) = C_renal_CV_timescale *C_rsna*(renal_sympathetic_nerve_activity - rsna_delayed);

###Disease effects (turned off by default)
#Glomerular hypertrophy
d/dt(disease_effects_increasing_Kf) = GP_effect_increasing_Kf;
#Loss of CD pressure natriuresis response over time
d/dt(disease_effects_decreasing_CD_PN) = CD_PN_loss_rate;

#Tubular hypertrophy
d/dt(tubular_length_increase) = PT_Na_reabs_effect_increasing_tubular_length;
d/dt(tubular_diameter_increase) = PT_Na_reabs_effect_increasing_tubular_diameter;


d/dt(water_out_s1_delayed) = C_renal_CV_timescale * C_pt_water*(water_out_s1 - water_out_s1_delayed);
d/dt(water_out_s2_delayed) = C_renal_CV_timescale * C_pt_water*(water_out_s2 - water_out_s2_delayed);
d/dt(water_out_s3_delayed) = C_renal_CV_timescale * C_pt_water*(water_out_s3 - water_out_s3_delayed);
d/dt(reabsorbed_urea_cd_delayed) = 0;#C_pt_water*(reabsorbed_urea_cd - reabsorbed_urea_cd_delayed);

#Urinary glucose excretion
d/dt(UGE) = C_renal_CV_timescale * RUGE;

#Serum Creatinine
d/dt(serum_creatinine) = C_renal_CV_timescale*(creatinine_synthesis_rate - creatinine_clearance_rate);
d/dt(cumNaExcretion) = C_renal_CV_timescale*Na_excretion_via_urine;
d/dt(cumWaterExcretion) = C_renal_CV_timescale*urine_flow_rate;
d/dt(cumCreatinineExcretion) = C_renal_CV_timescale*creatinine_clearance_rate;
d/dt(RTg_compensation) = C_renal_CV_timescale*excess_glucose_increasing_RTg;
d/dt(SGLT2_inhibition_delayed) = C_renal_CV_timescale*C_sglt2_delay*(SGLT2_inhibition - SGLT2_inhibition_delayed);
d/dt(RUGE_delayed) = C_renal_CV_timescale*C_ruge*(RUGE - RUGE_delayed);
d/dt(postglomerular_pressure_delayed) = C_renal_CV_timescale*C_postglomerular_pressure*(postglomerular_pressure - postglomerular_pressure_delayed); #necessary to prevent large oscillations
d/dt(postglomerular_pressure_error) = C_renal_CV_timescale*(postglomerular_pressure - RIHP0);
d/dt(renal_flow_rate_error) = C_renal_CV_timescale*(renal_blood_flow_L_min - nom_renal_blood_flow_L_min);
d/dt(MAP_delayed) = C_renal_CV_timescale*C_cycle2*(mean_arterial_pressure_MAP - MAP_delayed); #necessary to prevent large oscillations;
d/dt(RIHP_delayed)=C_renal_CV_timescale*C_cycle2*(RIHP - RIHP_delayed);
d/dt(Net_oncotic_pressure_diff) = C_renal_CV_timescale*C_cycle2*(Net_oncotic_pressure - Net_oncotic_pressure_diff);
d/dt(RISF) = tubular_reabsorption - capillary_filtration;
")

    test_that("large models compile", {
      expect_true(inherits(mod, "RxODE"))
    })
  },
  test = "lvl2"
)