Replies: 5 comments 14 replies
-
Hello jfndiaye Can you attach a script with the code you are trying to run? It will make it a lot easier to debug. |
Beta Was this translation helpful? Give feedback.
-
Codes are in the .txt files. |
Beta Was this translation helpful? Give feedback.
-
@jfndiaye, Do you know if the model ran on |
Beta Was this translation helpful? Give feedback.
-
Hi @jfndiaye I also am having an issue reproducing the simulation with the current It may take some time to figure out the differences between the version of One way to figure out what version that is is to check the R version that the paper is using (R 3.1.2) which predates the CRAN release of RxODE. So using the first released version of RxODE by: devtools::install_version("RxODE", version = "0.5-1") With this version the script works. Expand the details to see the code and results #library(rxode2)
#devtools::install_version("RxODE", version = "0.5-1")
library(RxODE)
library(minpack.lm) # library for least squares fit using levenberg-marquart algorithm
### From model_struct.R
############################## Define model structure ##########################
ode <- "
#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);
######################Systemic Hemodynamics #####################################
######Calculation of Systemic Vascular Resistance
#Systemic vascular resistance is a nominal value modulated by AngII and by a regulatory signal for tissue autoregulation necessary to maintain constant organ blood flow
###Whole body autoregulation mechanism wherein TPR adjusts to maintain constant organ blood flow (and thus constant cardiac output)
#Modeled as Proportional-Integral controller of TPR, where the input signal is the cardiac output error signal
tissue_autoregulation_signal = max(0.1,1+tissue_autoreg_scale*((Kp_CO/CO_scale_species)*(cardiac_output_delayed - CO_nom)+(Ki_CO/CO_scale_species)*CO_error));
###Effect of the RAAS (AT1-bound AngII) on systemic vascular resistance. For now, the slope is set to zero, i.e. AngII does not directly affect SVR
AT1_svr_int = 1 - AT1_svr_slope*nominal_equilibrium_AT1_bound_AngII;
AT1_bound_AngII_effect_on_SVR = AT1_svr_int + AT1_svr_slope * AT1_bound_AngII;
###Effect of RSNA on systemic vascular resistance
rsna_svr_int = 1 - rsna_svr_slope;
rsna_effect_on_svr = rsna_svr_int + rsna_svr_slope*rsna_delayed;
systemic_arterial_resistance = nom_systemic_arterial_resistance*tissue_autoregulation_signal*AT1_bound_AngII_effect_on_SVR*rsna_effect_on_svr;
###### Determination of Cardiac Output and Mean Arterial Pressure
#Cardiac output is a function of blood volume and resistance to venous return
###Empirical relationship between SVR, venous resistance, and resistance to venous return
resistance_to_venous_return = ((8 * R_venous + systemic_arterial_resistance) / 31);
###Empirical relationship between blood volume and mean filling pressure, taken from Karaaslan 2005; must be scaled for different species
mean_filling_pressure = nom_mean_filling_pressure + (blood_volume_L/BV_scale_species-blood_volume_nom)/venous_compliance;
venous_return = ((mean_filling_pressure) / resistance_to_venous_return);
cardiac_output = venous_return;
total_peripheral_resistance = systemic_arterial_resistance + R_venous;
mean_arterial_pressure_MAP = (cardiac_output * total_peripheral_resistance);
###Empirical relationship between blood volume and right atrial pressure, taken from Karaaslan 2005; must be scaled for different species
#Right atrial pressure is used as a signal for ANP and RSNA. For this purpose, this relationship is fine. However, a more complex model would be needed to truly model right atrial pressure
right_atrial_pressure = (nom_right_atrial_pressure * exp(blood_volume_L - 5*BV_scale_species));
###### Calculation of RSNA, ANP, and MR-bound Aldosterone
###Renal sympathetic nerve activity is assumed to be driven by changes in MAP and Right atrial pressure. The dominant effect is MAP, and the effect of RAP is much less, at least in the normal physiologic range.
MAP_effect_on_rsna = 0.5 + 1.0 / (1 + exp((mean_arterial_pressure_MAP - nominal_map_setpoint) / map_rsna_slope));
R_atrial_pressure_effect_on_rsna_int = 1+rap_rsna_slope*nom_right_atrial_pressure;
R_atrial_pressure_effect_on_rsna = R_atrial_pressure_effect_on_rsna_int - rap_rsna_slope * right_atrial_pressure;
renal_sympathetic_nerve_activity = nom_rsna*MAP_effect_on_rsna * R_atrial_pressure_effect_on_rsna*BB_effect_on_RSNA;
###ANP release is driven by changes in right atrial pressure
rap_anp_intercept = 1 + rap_anp_slope/2;
normalized_atrial_NP_concentration = nom_ANP*((rap_anp_intercept - rap_anp_slope / (1 + exp(right_atrial_pressure - nom_right_atrial_pressure))));
###Aldosterone bound to the mineralocorticoid receptor - normalized level. Aldosterone secretion is described later on.
aldosterone_concentration = normalized_aldosterone_level* nominal_aldosterone_concentration;
Aldo_MR_normalised_effect = normalized_aldosterone_level*MR_antagonist_effect_on_aldo_MR;
###################### Renal Vasculature #####################################
###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));
### 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
#The resistance of the arcuate, interlobular arterioles, and other vasculature prior the afferent arterioles is represented by a single resistance - the preafferent arteriole resistance
#The preafferent arterioles respond myogenically to changes in pressure, and also responds to AT1-bound AngII, RSNA, and ANP
#The dilation/constriction of the arterioles is limited, and thus the total combined effect of all regulators must saturate
preaff_arteriole_signal_multiplier = AT1_effect_on_preaff*rsna_effect_on_preaff*ANP_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 = 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
nom_afferent_arteriole_resistance = L_m3*viscosity_length_constant/(nom_afferent_diameter^4);
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_adjusted_signal_multiplier = (1/(1+exp(afferent_signal_nonlin_scale*(1-afferent_arteriole_signal_multiplier)))+0.5);
afferent_arteriole_resistance = nom_afferent_arteriole_resistance*afferent_arteriole_adjusted_signal_multiplier;
###### 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
nom_efferent_arteriole_resistance = L_m3*viscosity_length_constant/(nom_efferent_diameter^4);
efferent_arteriole_signal_multiplier = AT1_effect_on_eff * ANP_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 = 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;
temp=glomerular_pressure/(nom_glomerular_pressure+2);
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;
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;
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;
ECF_Na_concentration = ECF_sodium_amount/extracellular_fluid_volume;
sodium_storate_rate = Q_Na_store*((max_stored_sodium - stored_sodium)/max_stored_sodium)*(ECF_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_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((postglomerular_pressure- RIHP0) / pressure_natriuresis_PT_slope)));
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 - 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 - 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((postglomerular_pressure - RIHP0) / pressure_natriuresis_CD_slope)));
###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));
### 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 = rsna_PT_int + rsna_PT_scale/(1+exp((1 - rsna_delayed2)/rsna_PT_slope));
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
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 * 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;
#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(1,nominal_cd_na_reabsorption*cd_multiplier);
##################################### 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;
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_SGLT2_Na_reabs = SGTL2_Na_reabs_mmol_s1+SGTL2_Na_reabs_mmol_s2+SGTL1_Na_reabs_mmol; #mmol/min
#Under normal conditions, 3.9% of Na absorbed through SGLT2 (total_SGLT2_Na_reabs/SN_filtered_Na_load)
#Thus, this amount is subtracted to obtain the adjusted fractional reabsorption rate separate from SGLT2
e_pt_sodreab_adj = e_pt_sodreab*osmotic_natriuresis_effect_pt - 0.039;
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_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*exp(-Na_reabs_per_unit_length*L_pt_s1) - 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;
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;
##################################### 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;
#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;
################################################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;
ADH_water_permeability_old = min(1,max(0,nom_ADH_water_permeability*normalized_vasopressin_concentration));
ADH_water_permeability = normalized_vasopressin_concentration/(0.15+normalized_vasopressin_concentration);
#Water reabsorption follows gradient but is regulated by ADH
#max_water_reabs_CD = water_in_CD-(osmolality_in_CD*water_in_CD-2*Na_reabsorbed_CD)/osmolality_in_AscLoH;
#water_out_CD_before_osmotic_reabsorption = max(0,water_in_CD - ADH_water_permeability*max_water_reabs_CD);
osmoles_out_CD = osmoles_in_CD-2*(Na_in_CD - Na_out_CD);
#osmolality_out_CD_before_osmotic_reabsorption = osmoles_out_CD/water_out_CD_before_osmotic_reabsorption;
#water_reabsorbed_CD = ADH_water_permeability*water_out_CD_before_osmotic_reabsorption*(1-osmolality_out_CD_before_osmotic_reabsorption/osmolality_out_DescLoH);
#water_out_CD = water_out_CD_before_osmotic_reabsorption-water_reabsorbed_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;
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;
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;
#The amount of fluid reabsorbed is the difference between the amount filtered and the amount excreted
tubular_reabsorption = GFR_ml_min/1000 - urine_flow_rate;
#Renal Interstitial Hydrostatic Pressure
RIHP = postglomerular_pressure - (oncotic_pressure_peritubular_avg - interstitial_oncotic_pressure) + 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 = 4.9*mmHg_Nperm2_conv;#(renal_interstitial_hydrostatic_pressure-5)*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);
# Drug effects on RAAS pathways
ACEI_effect_on_ACE_activity = 1;
#
ARB_effect_on_AT1_binding = 1;
DRI_effect_on_PRA = 1;
###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
#e.g. Isaksson 2014
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));
plasma_renin_activity = concentration_to_renin_activity_conversion_plasma* plasma_renin_concentration*DRI_effect_on_PRA;
renin_secretion_rate = (log(2)/renin_half_life)*nominal_equilibrium_PRC*AT1_bound_AngII_effect_on_PRA*md_effect_on_renin_secretion*rnsa_effect_on_renin_secretion*HCTZ_effect_on_renin_secretion*renin_hyperactivity;
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;
ACE_activity = nominal_ACE_activity*ACEI_effect_on_ACE_activity;
chymase_activity = nominal_chymase_activity;
#AT1_receptor_binding_rate =nominal_AT1_receptor_binding_rate;
AT1_receptor_binding_rate = nominal_AT1_receptor_binding_rate*ARB_effect_on_AT1_binding;
AT2_receptor_binding_rate = nominal_AT2_receptor_binding_rate;
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 Extracellular fluid volume over time is determined by the different between water intake and urine outflow
d/dt(blood_volume_L) = C_water_intake_ecf_volume * (water_intake) + C_urine_flow_ecf_volume * (urine_flow_rate) + Q_water*(Na_concentration - ECF_Na_concentration);
d/dt(extracellular_fluid_volume) = Q_water*(ECF_Na_concentration - Na_concentration);
#Change in total body sodium over time is determined by the different between sodium intake and excretion
d/dt(sodium_amount) = C_na_excretion_na_amount * (Na_excretion_via_urine) + C_na_intake_na_amount * (Na_intake_rate) + Q_Na*(ECF_Na_concentration - Na_concentration);
d/dt(ECF_sodium_amount) = Q_Na*(Na_concentration - ECF_Na_concentration) - sodium_storate_rate;
d/dt(stored_sodium) = 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_tgf * (tubulo_glomerular_feedback_signal-tubulo_glomerular_feedback_effect);
d/dt(normalized_aldosterone_level) = C_aldo_secretion * (N_als-normalized_aldosterone_level);
d/dt(preafferent_pressure_autoreg_signal) = 500*(preafferent_pressure_autoreg_function - preafferent_pressure_autoreg_signal);
d/dt(glomerular_pressure_autoreg_signal) = 500*(glomerular_pressure_autoreg_function - glomerular_pressure_autoreg_signal);
d/dt(F_out_dt_delay) = 0;#100*(F_out_dt - F_out_dt_delay);
d/dt(cardiac_output_delayed) = C_cardiac_output_delayed*(cardiac_output - cardiac_output_delayed);
#Error signals for PI controllers of cardiac output and sodium concentration
d/dt(CO_error) = C_co_error*(cardiac_output-CO_nom);
d/dt(Na_concentration_error) = 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_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_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_P_bowmans*(P_in_pt_s1_mmHg - P_bowmans);
d/dt(oncotic_pressure_difference) = C_P_oncotic*(oncotic_pressure_avg - oncotic_pressure_difference);
d/dt(renal_blood_flow_L_min_delayed)=C_rbf*(renal_blood_flow_L_min - renal_blood_flow_L_min_delayed);
d/dt(renal_interstitial_hydrostatic_pressure) = C_rihp*(RIHP - renal_interstitial_hydrostatic_pressure);
d/dt(SN_macula_densa_Na_flow_delayed) = C_md_flow*( SN_macula_densa_Na_flow - SN_macula_densa_Na_flow_delayed);
d/dt(rsna_delayed) = C_rsna*(renal_sympathetic_nerve_activity - rsna_delayed);
d/dt(rsna_delayed2) = C_rsna2*(renal_sympathetic_nerve_activity - rsna_delayed2);
###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_pt_water*(water_out_s1 - water_out_s1_delayed);
d/dt(water_out_s2_delayed) = C_pt_water*(water_out_s2 - water_out_s2_delayed);
d/dt(water_out_s3_delayed) = 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) = RUGE;
#Serum Creatinine
d/dt(serum_creatinine) = creatinine_synthesis_rate - creatinine_clearance_rate;
d/dt(cumNaExcretion) = Na_excretion_via_urine;
d/dt(cumWaterExcretion) = urine_flow_rate;
d/dt(cumCreatinineExcretion) = creatinine_clearance_rate;
d/dt(RTg_compensation) = excess_glucose_increasing_RTg;
d/dt(SGLT2_inhibition_delayed) = C_sglt2_delay*(SGLT2_inhibition - SGLT2_inhibition_delayed);
d/dt(RUGE_delayed) = C_ruge*(RUGE - RUGE_delayed);
"
mod1 <- RxODE(model = ode, modName = "model1")
#load basecase
calcNomParams_human <- function(){
#Constants and Unit conversions
nL_mL=1e+06
dl_ml=0.01
L_dL=10
L_mL=1000
L_m3=0.001
g_mg=0.001
ng_mg=1e-06
secs_mins=60
min_hr=60
hr_day=24
min_day=1440
MW_creatinine=113.12
Pi=3.1416
pi=3.14
viscosity_length_constant=1.5e-09
gamma = 1.16667e-5; # viscosity of tubular fluid
mmHg_Nperm2_conv = 133.32
glucose_mg_mmol = 0.0056
uric_acid_mg_mmol = 0.006
#Scaling parameters - can be used to parameterize model for other species
ECF_scale_species = 1
BV_scale_species=1
water_intake_species_scale = 1
CO_scale_species = 1
########################################################################################
#Parameters of normal human physiology based on literature and commmon medical knowledge
########################################################################################
####Systemic parameters
nominal_map_setpoint=85 #mmHg
CO_nom= 5 #L/min
ECF_nom = 15
blood_volume_nom = 5 #L
Na_intake_rate=.07 #mEq/min - 100mmol/day or 2300 mg/day
nom_water_intake = 2.1 #L/day
ref_Na_concentration=138 #mEq/L
plasma_protein_concentration = 7 #g/dl
plasma_albumin_concentration= 35 #mg/ml
glucose_concentration = 5.5 #mmol/L
plasma_urea = 5 #mmol/L
nom_serum_uric_acid_concentration = 5.5 #mg/dl
equilibrium_serum_creatinine=0.92 #mg/dl
potassium_concentration=5 #mEq/L
P_venous=4 #mmHg
R_venous=3.4 #mmHg
nom_right_atrial_pressure=0.87 #mmHg
nom_mean_filling_pressure=7 #mmHg
venous_compliance = 1.2
####Renal parameters
nom_renal_blood_flow_L_min = 1 #L/min
baseline_nephrons=2e6
nom_Kf=3.9 #nl/min*mmHg
nom_oncotic_pressure_difference = 28.15 #mmHg
P_renal_vein=4 #mmHg
nom_oncotic_pressure_peritubular= 28.15 #mmHg
interstitial_oncotic_pressure = 5 #mmHg
#Renal Vasculature
nom_preafferent_arteriole_resistance= 14 ##15 #mmHg
nom_afferent_diameter=1.65e-5 ###1.5e-05 #mmHg
nom_efferent_diameter=1.1e-05 #mmHg
#Renal Tubules
Dc_pt_nom = 27e-6 #m
Dc_lh = 17e-6 #m
Dc_dt = 17e-6 #m
Dc_cd = 17e-6
L_pt_s1_nom = 0.005
L_pt_s2_nom = 0.005 #m
L_pt_s3_nom =0.004 #m
L_lh_des = 0.01 ###0.003 #m
L_lh_asc = 0.01 ###0.003 #m
L_dct = 0.005 #m
L_cd = L_lh_des
tubular_compliance = 0.2
Pc_pt_s1_mmHg = 20.2#19.4#13.2 #15 #mmHg
Pc_pt_s2_mmHg = 15
Pc_pt_s3_mmHg = 11 #mmHg
Pc_lh_des_mmHg = 8 #mmHg
Pc_lh_asc_mmHg = 7 #mmHg
Pc_dt_mmHg = 6 #mmHg
Pc_cd_mmHg = 5 #mmHg
nominal_pt_na_reabsorption=0.7# 0.75 #fraction
nominal_loh_na_reabsorption = 0.8 #0.65 #fraction
nominal_dt_na_reabsorption=0.5 #fraction
LoH_flow_dependence = 0.75
###Renal Glucose reabsorption
nom_glucose_reabs_per_unit_length_s1 = 5.4e-5
nom_glucose_reabs_per_unit_length_s2 = 0
nom_glucose_reabs_per_unit_length_s3 = 2.8e-5
diabetic_adaptation = 1
maximal_RTg_increase = 0.3
T_glucose_RTg = 100000
Anhe3 = 0
glucose_natriuresis_effect_pt = 0
glucose_natriuresis_effect_cd = 0
glucose_diuresis_effect_pt = 0
glucose_diuresis_effect_cd = 0
###Renal urea reabsorption
urea_permeability_PT = 0.5
####Renal albumin seiving
nom_glomerular_albumin_sieving_coefficient = 0.0001 #%
max_albumin_reabsorption_fraction=0.98448 #%
maximum_reabsorption_capacity = 5.5e-7 #mg/min
####RAAS Pathway parameters
concentration_to_renin_activity_conversion_plasma = 61
nominal_equilibrium_PRA = 1000 #937 #fmol/ml/hr
nominal_equilibrium_AngI = 7.5 #fmol/ml
nominal_equilibrium_AngII = 4.75 #fmol/ml
nominal_renin_half_life = 0.1733 # (hr)
nominal_AngI_half_life = 0.5/60 #(hr)
nominal_AngII_half_life = 0.66/60 #(hr)
nominal_AT1_bound_AngII_half_life = 12/60 #hr
nominal_AT2_bound_AngII_half_life = 12/60 #hr
ACE_chymase_fraction = 0.95 #% of AngI converted by ACE. The rest is converted by chymase
fraction_AT1_bound_AngII = 0.75 #assume AngII preferentially binds to AT1 vs AT2
########################################################################################
#The following parameters are calculated at equilibrium using the parameters above
########################################################################################
#This pressure is the setpoint that determines the myogenic response of the preafferent vasculature
nom_preafferent_pressure = nominal_map_setpoint - nom_renal_blood_flow_L_min*nom_preafferent_arteriole_resistance;
#This pressure is the setpoint that determines the myogenic response of the afferent vasculature
nom_glomerular_pressure = nom_preafferent_pressure - nom_renal_blood_flow_L_min*(L_m3*viscosity_length_constant/(nom_afferent_diameter^4)/baseline_nephrons);
#This pressure is the setpoint that determines the tubular pressure-natriuresis response
nom_postglomerular_pressure = nom_preafferent_pressure - nom_renal_blood_flow_L_min*(L_m3*viscosity_length_constant*(1/(nom_afferent_diameter^4)+1/(nom_efferent_diameter^4))/baseline_nephrons);
RIHP0 = nom_postglomerular_pressure
# The rate of sodium excretion must equal the rate of sodium intake. Sodium reabsorption rates vary along the tubule, but based on literature
# measurements we have a good, and literature data provides estimates for these rates. However, there is a precise
# rate of sodium reabsorption required to achieve the equilibrium defined by the parameters above.
# Assuming that reabsorption rates are known in all but one segment of the tubule, the exact rate
# of reabsorption of the remaining segment can be calculated. We chose to calculate the CD rate of reabsorpion based on estimates for
# PT, LoH, and DT reabsorption.
nom_GFR = nom_Kf*(nom_glomerular_pressure - nom_oncotic_pressure_difference - Pc_pt_s1_mmHg)/nL_mL*baseline_nephrons;
nom_filtered_sodium_load = nom_GFR/L_mL*ref_Na_concentration;
nom_PT_Na_outflow = nom_filtered_sodium_load*(1-nominal_pt_na_reabsorption);
nom_Na_in_AscLoH = nom_PT_Na_outflow/baseline_nephrons;
AscLoH_Reab_Rate =(2*nominal_loh_na_reabsorption*nom_Na_in_AscLoH)/L_lh_des; #osmoles reabsorbed per unit length per minute. factor of 2 because osmoles = 2
nom_LoH_Na_outflow = nom_PT_Na_outflow*(1-nominal_loh_na_reabsorption);
nom_DT_Na_outflow = nom_LoH_Na_outflow*(1-nominal_dt_na_reabsorption);
nominal_cd_na_reabsorption = 1-Na_intake_rate/nom_DT_Na_outflow;
#RBF = (MAP - P_venous)/RVR. Given MAP, P_venous, RBF, and preafferent, afferent, and efferent resistances, the remaining peritubular resistance at steady state can be determined
nom_RVR = (nominal_map_setpoint - P_venous)/nom_renal_blood_flow_L_min
nom_peritubular_resistance = nom_RVR - (nom_preafferent_arteriole_resistance + L_m3*viscosity_length_constant*(1/nom_afferent_diameter^4+1/nom_efferent_diameter^4)/baseline_nephrons);
#Calculate the normal amount of sodium reabsorbed per unit surface area of the PT
PT_Na_reab_perUnitSA_0 = (nom_filtered_sodium_load/baseline_nephrons)* nominal_pt_na_reabsorption/(3.14*Dc_pt_nom*(L_pt_s1_nom+L_pt_s2_nom+L_pt_s3_nom))
#Given the values for baseline MAP and CO above, the baseline TPR required to maintain this MAP and CO can be calculated. Since TPR includes renal vascular resistance, the baseline systemic (non-renal) resistance
#can be calculated from this TPR and the values for baseline renal resistances defined above.
nom_TPR = nominal_map_setpoint/CO_nom
#nom_systemic_arterial_resistance= (nom_TPR-R_venous)*nom_RVR/(nom_RVR - nom_TPR-R_venous)
nom_systemic_arterial_resistance= nom_TPR-R_venous
#Calculation of peritubular ultrafiltration coefficient
tubular_reabsorption = nom_GFR/1000 - nom_water_intake*water_intake_species_scale/60/24 #at SS, water excretion equals water intake
#Both RIHP and Kf are unknown, so we can either assume RIHP and calculate Kf, or vice versa. Since RIHP has been measured experimentally,
#it seems better to assume a normal value for RIHP and calculate Kf
nom_peritubular_cap_Kf = - tubular_reabsorption/(nom_postglomerular_pressure - RIHP0 - (nom_oncotic_pressure_peritubular - interstitial_oncotic_pressure))
#Creatinine synthesisrate at equilibrium
creatinine_synthesis_rate = equilibrium_serum_creatinine * dl_ml * nom_GFR #Units: mg/min
#### Uric Acid reabsorption
uric_acid_fractional_reabs_rate = 0.95 #%
uric_acid_synthesis_rate = nom_serum_uric_acid_concentration * dl_ml * nom_GFR*(1-uric_acid_fractional_reabs_rate) #Units: mg/min
####RAAS Pathway parameters
#Values for half lives and equilibrium concentrations of RAAS peptides available in the literature and
# defined above to calculate nominal values for other RAAS parameters not available in the literature:
#ACE activity
#Chymase activity
#AT1 receptor binding rate
#AT2 receptor binding rate
#equilibrium AT1_bound_AngII
#These values are then assumed to be fixed unless specified otherwise.
#Calculating these nominal parameter values initially in a separate file is required so that these parameters can then be varied independently in the main model
nominal_equilibrium_PRC = nominal_equilibrium_PRA/concentration_to_renin_activity_conversion_plasma
nominal_AngI_degradation_rate = log(2)/nominal_AngI_half_life #/hr
nominal_AngII_degradation_rate = log(2)/nominal_AngII_half_life #/hr
nominal_AT1_bound_AngII_degradation_rate = log(2)/nominal_AT1_bound_AngII_half_life
nominal_AT2_bound_AngII_degradation_rate = log(2)/nominal_AT2_bound_AngII_half_life
#ACE converts 95% of AngI, chymase converts the rest
nominal_ACE_activity = (ACE_chymase_fraction*(nominal_equilibrium_PRA - nominal_AngI_degradation_rate*nominal_equilibrium_AngI)/nominal_equilibrium_AngI)#Therapy_effect_on_ACE
nominal_chymase_activity = (1-ACE_chymase_fraction)*(nominal_equilibrium_PRA - nominal_AngI_degradation_rate*nominal_equilibrium_AngI)/nominal_equilibrium_AngI
#75% of bound AngII is AT1, the rest is AT2
nominal_AT1_receptor_binding_rate = fraction_AT1_bound_AngII*(nominal_equilibrium_AngI*(nominal_ACE_activity+nominal_chymase_activity)-nominal_AngII_degradation_rate*nominal_equilibrium_AngII)/nominal_equilibrium_AngII
nominal_AT2_receptor_binding_rate = (1-fraction_AT1_bound_AngII)*(nominal_equilibrium_AngI*(nominal_ACE_activity+nominal_chymase_activity)-nominal_AngII_degradation_rate*nominal_equilibrium_AngII)/nominal_equilibrium_AngII
nominal_equilibrium_AT1_bound_AngII = nominal_equilibrium_AngII*nominal_AT1_receptor_binding_rate/nominal_AT1_bound_AngII_degradation_rate
nominal_equilibrium_AT2_bound_AngII = nominal_equilibrium_AngII*nominal_AT2_receptor_binding_rate/nominal_AT2_bound_AngII_degradation_rate
########################################################################################
#The following parameters were determined indirectly from many different literature studies on the response
#various changes in the system (e.g. drug treatments, infusions of peptide, fluid, sodium, etc.....)
########################################################################################
#Effects of AT1-bound AngII on preafferent, afferent, and efferent resistance, and aldosterone secretion
AT1_svr_slope = 0
AT1_preaff_scale = 1 #0.4
AT1_preaff_slope = 7
AT1_aff_scale=0.5
AT1_aff_slope=7
AT1_eff_scale=0.1
AT1_eff_slope=7
AT1_PT_scale = 0.075
AT1_PT_slope = 7
AT1_aldo_slope = 0.01
#Effects of Aldosterone on distal and collecting duct sodium reabsorption
nominal_aldosterone_concentration=85
K_Na_ratio_effect_on_aldo = 1
aldo_DCT_scale=0.08
aldo_DCT_slope = 0.5
aldo_CD_scale=0.35
aldo_CD_slope = 0.5
#Effects of Atrial Natriuretic Peptide (ANP)preafferent, afferent, and efferent resistance and collectin duct sodium reabsorption
nom_ANP=1
rap_anp_slope=1
ANP_preaff_scale = 2
ANP_preaff_slope = 1
ANP_aff_scale = 2
ANP_aff_slope = 1
ANP_eff_scale = .4
ANP_eff_slope = 1
anp_CD_scale =0.25
anp_CD_slope = 1
#Effects of Renal Sympathetic Nerve Activity on preafferent resistance, renin secretion, and PT sodium reabsorption
nom_rsna = 1
map_rsna_slope=5
rap_rsna_slope=0.008
rsna_preaff_scale = 0
rsna_preaff_slope = 0.25
rsna_PT_scale=0
rsna_PT_slope=1
rsna_CD_scale = 0
rsna_CD_slope = 1
rsna_renin_slope=0.36
rsna_svr_slope = 0
#Na and water transfer between blood, ECF
Q_water = 1
Q_Na = 1
Q_Na_store = 0 #80
max_stored_sodium = 500 #meq
#Osmolarity control of vasopressin secretion
Na_controller_gain=0.05
Kp_VP = 1
Ki_VP = 0.01
nom_ADH_urea_permeability = .98
nom_ADH_water_permeability = .98
#Effects of Vasopressin on water intake and reabsorption
nominal_vasopressin_conc=4
water_intake_vasopressin_scale = 0.25#1.5
water_intake_vasopressin_slope = -0.5
#Magnitude and Steepness of tubuloglomerular feedback
S_tubulo_glomerular_feedback=0.7
F_md_scale_tubulo_glomerular_feedback=6
MD_Na_concentration_setpoint = 59.7
#Effect of macula densa sodium flow on renin secretion
#md_renin_low_slope = 5
#md_renin_high_slope = -0.2
md_renin_A = 1
md_renin_tau = 1.25
renin_hyperactivity = 1
#Responsiveness of renal vasculature to regulatory signals
preaff_diameter_range=0.25
afferent_diameter_range=1.2e-05
efferent_diameter_range=3e-06
preaff_signal_nonlin_scale=3
afferent_signal_nonlin_scale=3
efferent_signal_nonlin_scale=3
#Limit on PT sodium reabsorption
renal_threshold_Na_reabs = 16e-6
#Empirical relationship between blood volume and cardiac filling pressure - from Guyton
BV_filling_pressure_slope=7.436
#RAAS pathway (these parameters can be set to different values than used to calculate the equilibrium state above)
AngI_half_life=0.008333
AngII_half_life=0.011
AT1_bound_AngII_half_life=0.2
AT1_PRC_slope=-0.9 #1.2
AT1_PRC_yint=0
AT2_bound_AngII_half_life=0.2
concentration_to_renin_activity_conversion_plasma=61
fraction_AT1_bound_AngII=0.75
nominal_ACE_activity=48.9
nominal_AT1_receptor_binding_rate=12.1
nominal_AT2_receptor_binding_rate=4.0
nominal_chymase_activity=1.25
nominal_equilibrium_AT1_bound_AngII=16.6315288606173
nominal_equilibrium_PRC=16.4
renin_half_life=0.1733
#Transfer constants for ODEs - determine speed of processes
C_aldo_secretion=1000
C_prerenal_blood_pressure=1000
C_P_bowmans = 1000
C_P_oncotic = 1000
C_rbf=1000
C_pt_water=1000
C_rsna = 1000
C_rsna2 = 1
C_tgf_reset=0
C_cardiac_output_delayed=.001
C_co_error=0.00001
C_vasopressin_delay = 0.01
C_ruge = 0.0001
C_md_flow = 0.001#Time delay between MD sodium flow and renin secretion
C_rihp = 0.01#0.01 #time delay between peritubular pressure and RIHP
C_tgf=1#/30 #1000
C_na_excretion_na_amount=-1#/2#/60
C_na_intake_na_amount=1#/2#/60
C_urine_flow_ecf_volume=-1#/2#/60
C_water_intake_ecf_volume=1#/2#/60
C_Na_error=0.1 #1/6#0
C_serum_creatinine = 1#/2#/60
#Therapy effects
HCTZ_effect_on_DT_Na_reabs = 1
HCTZ_effect_on_renin_secretion = 1
DRI_effect_on_PRA = 1
CCB_effect_on_preafferent_resistance = 1
CCB_effect_on_afferent_resistance = 1
CCB_effect_on_efferent_resistance = 1
MR_antagonist_effect_on_aldo_MR = 1
####################################################################
#These parameters are by default set to ensure strong autoregulation of cardiac output, RBF, glomerular pressure, and MAP
#However, reducing these parameters reduces the ability of the system to autoregulate, and is necessary for modeling the development of hypertension, etc.
####################################################################
#Metabololic tissue autoregulation of cardiac output
tissue_autoreg_scale=0.2
Kp_CO=1.5
Ki_CO=30
#Renal autoregulation of glomerular pressure
gp_autoreg_scale=0
preaff_autoreg_scale = 0
myogenic_steepness=2
#Renal autoregulation of renal blood flow
RBF_autoreg_scale = 0#3
RBF_autoreg_steepness=0.001
#Pressure natiuresis effect through collecting duct sodium reabsorption
#Parameters selected based on Isaksson 2014:
#For a 10X increase in salt intake:MAP increases by 5mmHg, Renin decreases by 45%
#GFR increases by 1.4ml/min
#Strong CD effect required to minimize BP rise
#PT effect + LoH effect required to produce renin response
#If PT effect is too big, GFR will decrease instead of increase.
#So LoH must make up for the rest
max_pt_reabs_rate = 0.995
pressure_natriuresis_PT_scale = 0.1
pressure_natriuresis_PT_slope = 1
pressure_natriuresis_LoH_scale = 0
pressure_natriuresis_LoH_slope = 1
pressure_natriuresis_DCT_scale = 0
pressure_natriuresis_DCT_slope = 1
max_cd_reabs_rate = 0.995
pressure_natriuresis_CD_scale = 0.1
pressure_natriuresis_CD_slope=1
#Glomerular pressure effect on glomerular hypertrophy
maximal_glom_surface_area_increase = 0#.5
T_glomerular_pressure_increases_Kf = 2000
#PT sodium reabsorption effects on tubular hypertrophy
maximal_tubule_length_increase = 0#.5
maximal_tubule_diameter_increase = 0#.25
T_PT_Na_reabs_PT_length = 1e10
T_PT_Na_reabs_PT_diameter = 1e10
####Proteinuria
nom_glomerular_albumin_sieving_coefficient = 6.2e-4
Emax_seiving = 6000
Gamma_seiving = 2
Km_seiving = 30
max_PT_albumin_reabsorption_rate = 0.1
nom_albumin_excretion_rate = 3.5e-9
#reabsorptive capacity based on observation of no proteinuria after UNX
SN_albumin_reabsorptive_capacity = 1.4e-6 #mg/min/tubule
max_PT_albumin_reabsorption_rate=0.2
Emax_seiving=4#100
Km_seiving=25
Gamma_seiving=3
nom_GP_seiving_damage=65
c_albumin = 0.0231 #min/nl, from Dean and Lazzara
seiving_inf = 4.25e-4 #from Dean and Lazarra, calculated for seiving coeff =0.00062 when SNGFR = 50 nl/min
maximal_glom_surface_area_increase=0.5
#Reduce Kf due to glomerulosclerosis
disease_effects_decreasing_Kf = 0
#Disease effects
disease_effect_on_nephrons = 0
max_s1_Na_reabs = 7.5e-6
max_s2_Na_reabs = 1#3e-6
max_s3_Na_reabs = 1#2e-6
max_deltaLoH_reabs=0.75e-6#2e-6
CD_Na_reabs_threshold = 2.5e-7
#Rate at which the tubular pressure natriuresis mechanism is lost in diabetes (should be zero or negative number)
CD_PN_loss_rate = 0
#Treatment Parameters
BB_effect_on_RSNA = 1
SGLT2_inhibition = 1
SGLT1_inhibition = 1
C_sglt2_delay = 0.1
CA_inhibitor = 1
ACEi_effect_on_ACE = 1
loop_diuretic_effect = 1
t=sort(ls())
param=sapply(t,names)
for (i in 1:length(t)){
param[i]=get(t[i])
}
param$param=NULL
return(param)
}
theta=unlist(calcNomParams_human())
###Initial conditions - do NOT change order!!!
#Order must match order in model file
#labels are not used by RxODE to match init to compartment
inits <- c( AngI=8.164, AngII=5.17,
AT1_bound_AngII=16.6, AT2_bound_AngII = 5.5, plasma_renin_concentration=17.845,
blood_volume_L = theta["blood_volume_nom"],
extracellular_fluid_volume=theta["ECF_nom"],
sodium_amount= as.numeric(theta["blood_volume_nom"])*as.numeric(theta["ref_Na_concentration"]),
ECF_sodium_amount= as.numeric(theta["ECF_nom"])*as.numeric(theta["ref_Na_concentration"]),
stored_sodium = 0, #relative number - actual value not known
tubulo_glomerular_feedback_effect=1,
normalized_aldosterone_level=1,
preafferent_pressure_autoreg_signal=1,
glomerular_pressure_autoreg_signal=1,
F_out_dt_delay=5e-12,
cardiac_output_asdelayed=theta["CO_nom"],
CO_error=0, Na_concentration_error = 0,
normalized_vasopressin_concentration_delayed = 1,
F0_TGF=theta["nom_LoH_Na_outflow"],
P_bowmans=theta["Pc_pt_s1_mmHg"],
oncotic_pressure_difference=theta["nom_oncotic_pressure_difference"],
renal_blood_flow_L_min_delayed=theta["nom_renal_blood_flow_L_min"],
renal_interstitial_hydrostatic_pressure = theta["nom_postglomerular_pressure"],
SN_macula_densa_Na_flow_delayed = as.numeric(theta["nom_LoH_Na_outflow"])/as.numeric(theta["baseline_nephrons"]),
rsna_delayed = 1,
rsna_delayed2 = 1,
disease_effects_increasing_Kf= 0, disease_effects_decreasing_CD_PN = 0,
tubular_length_increase=0, tubular_diameter_increase=0,
Water_out_S1_delayed=3e-8,
Water_out_S2_delayed=1.9e-8,
Water_out_S3_delayed=1.2e-8,
reabsorbed_urea_cd_delayed =0,#10e-8,
water_out_DescLoH_delayed = 8e-8,
osmolality_in_AscLoH_delayed = 826,
UGE = 0,
serum_creatinine = as.numeric(theta["equilibrium_serum_creatinine"])*as.numeric(theta["blood_volume_nom"]), cumNaExcretion = 0, cumWaterExcretion = 0, cumCreatinineExcretion = 0,
cumNaExcretion = 0,
cumWaterExcretion = 0,
cumCreatinineExcretion = 0,
RTg_compensation = 0,
SGLT2_inhibition_delayed = 1,
RUGE_delayed = 0
)
########Run to equilibrium
ev <- eventTable()
ev$add.sampling(seq(0,100,by=1))
par(mfrow=c(1,3))
x <- mod1$run(theta, ev, inits=inits, atol = 1.0e-8, rtol = 1.0e-6,
hmin=0, hmax=0, hini=0, safeZero=FALSE,
maxordn = 12L, maxords = 5L, method="lsoda")
matplot(x[,"mean_arterial_pressure_MAP"],type="l")
plot(x[,"extracellular_fluid_volume"])
plot(x[,"GFR_ml_min"]) Created on 2022-07-26 by the reprex package (v2.0.1) |
Beta Was this translation helpful? Give feedback.
-
If you are trying in your mac, use the M1 CRAN supported compilers: |
Beta Was this translation helpful? Give feedback.
-
Hi everyone,
I'm new in the R computation area. Now I'm trying to run a cardiorenal model already made by Hallow et al by using rxode2 instead of RxODE (since I had trouble installing RxODE on my M1 mac). When I try to run the program, I get the following error message :
Any idea for this kind of problem would be appreciated.
Thanks
Beta Was this translation helpful? Give feedback.
All reactions