PBPK Sensitivity Analysis
Welcome to this tutorial on sensitivity analysis in physiologically-based pharmacokinetic (PBPK) modeling. The essence of PBPK models lies in their ability to anticipate drug behavior in the human body, accounting for variations in physiological, biological, and molecular properties. However, the reliability and utility of such models rest upon the robustness of their parameters. That's where sensitivity analysis comes into play.
Sensitivity analysis is a crucial method in PBPK modeling that allows us to understand how changes in model parameters can influence the model output. In essence, it gives us insights into which parameters have a significant impact on the model and which ones are less critical. By performing sensitivity analysis, we can not only enhance the credibility of our model but also identify the parameters that require precise estimation or further investigation.
This tutorial is designed for both beginners and intermediate users who have a basic understanding of pharmacokinetics and PBPK modeling. We'll begin by introducing the concept of sensitivity analysis and why it's crucial in PBPK modeling. We'll then delve into various methods of performing sensitivity analysis, ranging from simple univariate approaches to more sophisticated global methods. We'll also discuss how to interpret the results from these analyses to refine your PBPK model and validate its predictive power.
Finally, we'll provide hands-on examples using common software tools in the field, and walk through the process of conducting a sensitivity analysis on a typical PBPK model. We believe that through this practical approach, you'll gain a solid grasp of this essential aspect of PBPK modeling, enabling you to make more informed decisions in your research or pharmaceutical development processes.
Stay tuned for a comprehensive journey through the intricacies of sensitivity analysis in PBPK modeling, and we hope you'll find this tutorial both informative and practical.
As we delve deeper into sensitivity analysis, we encounter various techniques for conducting these analyses. Among them, one technique particularly stands out for its efficiency and accuracy when dealing with high-dimensional models such as PBPK models - the global sensitivity analysis using Latin Hypercube Sampling (LHS).
- Global Sensitivity Analysis
Unlike local sensitivity analysis, which explores the impact of one parameter at a time while holding others constant, global sensitivity analysis investigates the influence of all parameters, over their entire range, considering their interactions. It provides a holistic view of how parameter uncertainties contribute to the model output variability. This global approach is critical in complex PBPK models where parameters often interact and influence each other, leading to non-linear, non-additive effects.
- Latin Hypercube Sampling (LHS)
One of the challenges in global sensitivity analysis is the selection of parameter sets that provide a thorough exploration of the parameter space. For complex models with a high number of parameters, the volume of the parameter space increases exponentially, rendering many sampling techniques inefficient.
Latin Hypercube Sampling (LHS) addresses this challenge. LHS is a stratified sampling method that ensures a more thorough and efficient exploration of the parameter space, reducing the computational burden compared to traditional methods. It operates by dividing the distribution of each parameter into equally probable intervals, or strata. A sample is then randomly drawn from each stratum, and the selected samples from each parameter are combined into sets in such a way that every parameter set contains exactly one randomly chosen value from each parameter. This strategy ensures that the entire range of each parameter is covered, leading to a well-distributed sample of the parameter space.
Applying LHS to global sensitivity analysis in PBPK models offers multiple benefits. It reduces the risk of overlooking important regions of the parameter space and provides a more robust estimate of how parameter uncertainties influence model output. Moreover, it helps in identifying which parameters drive the model behavior, thereby aiding in model refinement and optimization.
In the following sections of this tutorial, we'll explore how to implement global sensitivity analysis using Latin Hypercube Sampling in a PBPK model context, with hands-on examples and best practices for interpretation of results.
After carrying out a global sensitivity analysis using Latin Hypercube Sampling (LHS), the subsequent step involves interpreting the results. One of the powerful techniques for this purpose is the use of partial correlation coefficients.
- Understanding Partial Correlation Coefficients
Partial correlation is a measure of the relationship between two variables while controlling for the effect of one or more additional variables. In the context of PBPK modeling, a partial correlation coefficient measures the relationship between a model parameter (input variable) and the model output, while controlling for all other parameters.
- Interpreting Partial Correlation Coefficients in Sensitivity Analysis
Partial correlation coefficients can be used to rank the parameters according to their influence on the output. A high absolute value of the partial correlation coefficient for a parameter indicates that it has a substantial influence on the model output. If the coefficient is positive, it means that the parameter and the output are positively correlated (as the parameter increases, so does the output), and if it's negative, they are inversely correlated.
Furthermore, the sign of the coefficient can provide insights into the direction of the relationship. If the coefficient is positive, it implies that as the parameter value increases, the output increases. On the other hand, a negative coefficient signifies an inverse relationship, meaning that as the parameter value increases, the output decreases.
- Advantages and Limitations
Using partial correlation coefficients to interpret the results of a global sensitivity analysis can provide a robust way of understanding the interplay between various parameters and the output. However, it assumes a linear relationship between the parameters and the output. This assumption might not hold in complex PBPK models where parameters can interact in a nonlinear or non-monotonic fashion.
- Practical Application
In the next part of this tutorial, we will walk through a practical example of conducting a global sensitivity analysis in a PBPK model using Latin Hypercube Sampling. We will then illustrate how to calculate and interpret partial correlation coefficients to identify the key parameters driving model output variability.
Remember, the goal of sensitivity analysis is not merely to generate data but to enhance the model by understanding the implications of the results. Thus, careful interpretation is crucial to refining your PBPK model and improving its predictive performance.
References:
- Primary Source: Scotcher, D., Jones, C., Rostami-Hodjegan, A., & Galetin, A. (2016). Novel minimal physiologically-based model for the prediction of passive tubular reabsorption and renal excretion clearance. European Journal of Pharmaceutical Sciences, 94, 59–71. https://doi.org/10.1016/j.ejps.2016.03.018
- Secondary Source: Huang, W., & Isoherranen, N. (2018). Development of a Dynamic Physiologically Based Mechanistic Kidney Model to Predict Renal Clearance. CPT: Pharmacometrics & Systems Pharmacology, 7(9), 593–602. https://doi.org/10.1002/psp4.12321
- Whole-body PBPK: Tsamandouras, N., Kostrezewski, T., Stokes, C.L., Hughest, D.J., & Cirit, M. (2018). Quantitative Assessment of Population Variability in Hepatic Drug Metabolism Using a Perfused Three-Dimensional Human Liver Microphysiological System. Journal of Pharmacology and Experimental Therapeutics, 360, 95-105. http://dx.doi.org/10.1124/jpet.116.237495
% Define parameter search space for repeated runs and sensitivity analyses
num_samples = 100;
% convert parameters into vector of values to test
param_set_id = 'global';
switch param_set_id
case 'global'
% logP, fup, pKa, PappAB, CLintl
ionization = 'ACID'; % 'NEUTRAL','ACID','BASE'
varied = {'model_PkidAB','model_clint_hep','model_logP','model_fup','model_pKa'};
x = lhsdesign(num_samples,length(varied));
pmin = 0; % PkidAB
pmax = 50;
PkidAB = pmin + x(:,1)'*(pmax-pmin);
pmin = 0; % clint_hep
pmax = 150;
clint_hep = pmin + x(:,2)'*(pmax-pmin);
pmin = -5; % logP
pmax = 5;
logP = pmin + x(:,3)'*(pmax-pmin);
pmin = 0; % fup
pmax = 1;
fup = pmin + x(:,4)'*(pmax-pmin);
pmin = 2; % pKa
pmax = 12;
pKa = pmin + x(:,5)'*(pmax-pmin);
case 'PkidAB'
varied = {'model_PkidAB'};
x = lhsdesign(num_samples,length(varied));
pmin = 0; % x1e-6 cm/s
pmax = 30; % x1e-6 cm/s
PkidAB = pmin + x(:,1)'*(pmax-pmin); % cm^2
num_samples = length(PkidAB);
end
% Model equations
eqns = {
'model:'
% Drug-related parameters
% These need to be initialized in eqns placeholder values will not be
% used as long as they exist with the same name in the parent script,
% provided they are called in P = {...} at the end.
'PkidAB = 0' % 1x-6 cm/s, in vitro kidney permeability (apical->basal)
'PkidBA = 0' % 1x-6 cm/s, in vitro kidney permeability (basal-apical)
'fup = 1' % fraction unbound in the plasma
'Rbp = 1' %Blood-plasma partition coeff
'clint_hep = 0' %Hepatic intrinsic clearance
'fub = fup/Rbp' % fraction unbound in blood
%Initialize values needed for the Rogers and Rowland model.
'logP=1'
'pKa=1'
'Rbp=1'
'Ht=0.45'
'ionization=''NEUTRAL'''
%Extract partition coefficients for tissues using the Rogers and
%Rowland model. The 'rest of body' compartment is assumed to be an
%average of the other tissues.
'Kp_blood = tissueplasma_RR_blood_wrapper(ionization,pKa,logP,fup,Rbp,Ht)'
'Kpb_adipose = Kp_blood(1,:)'
'Kpb_bone = Kp_blood(2,:)'
'Kpb_brain = Kp_blood(3,:)'
'Kpb_gut = Kp_blood(4,:)'
'Kpb_heart = Kp_blood(5,:)'
'Kpb_kidney = Kp_blood(6,:)'
'Kpb_liver = Kp_blood(7,:)'
'Kpb_lung = Kp_blood(8,:)'
'Kpb_muscle = Kp_blood(9,:)'
'Kpb_skin = Kp_blood(10,:)'
'Kpb_spleen = Kp_blood(11,:)'
'Kpb_thymus = Kp_blood(12,:)'
'Kpb_rest = mean(Kp_blood(1:11,:),1)'
% Physiological parameters
'BW = 70' %body weight, in kg
'GFR = 120' %glomerular filtration rate, in ml/min
% Pull and extract physiological parameters associated with
% Tsamandouras et al, 2017.
'phys_params = get_physiological_parameters(BW)'
%Volumes (ml)
'vol_arterial = phys_params.volume.arterial'
'vol_lung = phys_params.volume.lung'
'vol_venous = phys_params.volume.venous'
'vol_spleen = phys_params.volume.spleen'
'vol_gut = phys_params.volume.gut'
'vol_liver = phys_params.volume.liver'
'vol_brain = phys_params.volume.brain'
'vol_heart = phys_params.volume.heart'
'vol_kidney = phys_params.volume.kidney'
'vol_skin = phys_params.volume.skin'
'vol_muscle = phys_params.volume.muscle'
'vol_adipose = phys_params.volume.adipose'
'vol_bone = phys_params.volume.bone'
'vol_rest = phys_params.volume.rest'
%Flow rates (ml/min)
'Q_arterial = phys_params.flow_rate.arterial'
'Q_lung = phys_params.flow_rate.lung'
'Q_venous = phys_params.flow_rate.venous'
'Q_spleen = phys_params.flow_rate.spleen'
'Q_gut = phys_params.flow_rate.gut'
'Q_liver = phys_params.flow_rate.liver'
'Q_brain = phys_params.flow_rate.brain'
'Q_heart = phys_params.flow_rate.heart'
'Q_kidney = phys_params.flow_rate.kidney'
'Q_skin = phys_params.flow_rate.skin'
'Q_muscle = phys_params.flow_rate.muscle'
'Q_adipose = phys_params.flow_rate.adipose'
'Q_bone = phys_params.flow_rate.bone'
'Q_rest = phys_params.flow_rate.rest'
% Pull and extract physiological parameters associated with
% Scotcher et al, 2016, for kidney regions
'scotcher_params = get_scotcher_parameters()'
%Kidney region pattern: (variable)_(vascular/tubular)_(region)
%Where regions are proximal tubule (pt), loop of Henle (loh),
%distal tubule (dt) or collecting duct (cd)
%Volumes (ml)
'vol_t_pt = scotcher_params.volume.proximal_tubule'
'vol_t_loh = scotcher_params.volume.loop_of_henle'
'vol_t_dt = scotcher_params.volume.distal_tubule'
'vol_t_cd = scotcher_params.volume.collecting_duct'
%Surface areas (cm^2)
'sa_t_pt = scotcher_params.surface_area.proximal_tubule'
'sa_t_loh = scotcher_params.volume.loop_of_henle'
'sa_t_dt = scotcher_params.volume.distal_tubule'
'sa_t_cd = scotcher_params.volume.collecting_duct'
%Flow rates (ml/min)
'Q_t_pt = scotcher_params.flow_rate.proximal_tubule'
'Q_t_loh = scotcher_params.flow_rate.loop_of_henle'
'Q_t_dt = scotcher_params.flow_rate.distal_tubule'
'Q_t_cd = scotcher_params.flow_rate.collecting_duct'
%Simplifying assumption: vascular kidney regions have identical volumes
%as their tubular counterparts (not sure how true this is for
%collecting ducts?)
'vol_v_pt = vol_t_pt'
'vol_v_loh = vol_t_loh'
'vol_v_dt = vol_t_dt'
'vol_v_cd = vol_t_cd'
% scaling factors (not all will necessarily be used)
'secpermin = 60' %s to min
'pappmulti = 10^(-6)' %1e-6 scaling for Papp
'papp_SF = secpermin.*pappmulti'
% Active transport
'Km = 1'
'Vmax = 0'
'MM(C) = Vmax*C/(C+Km);'
%Infusion equations
%Note the functional form in the rate equation essentially makes the
%bw_dosing_flag variable a trigger for multiplying by body weight; if
%it is zero, the dose is multiplied by 1, and if 1, multiplied by BW.
'infusion_duration = 10'
'bw_dosing_flag = 0' %Indicates whether dosing is per bw or not
'ton=0; toff=infusion_duration' %On and off for infusion time
'dose = 9999' %ng
'weight = BW' %kg
'rate = dose.*weight./toff'
'drug(t) = rate.*(t>=ton&t<=toff)'
%ODEs
%Transit compartments
'dC_lung/dt = (1/vol_lung).*(Q_lung.*(C_venous - C_lung./Kpb_lung)); C_lung(0) = 0*ones(1,Npop)' %Lung (receives input from venous compartment, outputs to arterial)
'dC_arterial/dt = (1/vol_arterial).*(Q_lung.*(C_lung./Kpb_lung - C_arterial)); C_arterial(0) = 0*ones(1,Npop)' %Arterial compartment (input from lung, outputs most other tissues)
'dC_spleen/dt = (1/vol_spleen).*(Q_spleen.*(C_arterial - C_spleen./Kpb_spleen)); C_spleen(0) = 0*ones(1,Npop)' %Spleen
'dC_gut/dt = (1/vol_gut).*(Q_gut.*(C_arterial - C_gut./Kpb_gut)); C_gut(0) = 0*ones(1,Npop)' %Gut
'dC_brain/dt = (1/vol_brain).*(Q_brain.*(C_arterial - C_brain./Kpb_brain)); C_brain(0) = 0*ones(1,Npop)' %Brain
'dC_heart/dt = (1/vol_heart).*(Q_heart.*(C_arterial - C_heart./Kpb_heart)); C_heart(0) = 0*ones(1,Npop)' %Heart
'dC_skin/dt = (1/vol_skin).*(Q_skin.*(C_arterial - C_skin./Kpb_skin)); C_skin(0) = 0*ones(1,Npop)' %Skin
'dC_muscle/dt = (1/vol_muscle).*(Q_muscle.*(C_arterial - C_muscle./Kpb_muscle)); C_muscle(0) = 0*ones(1,Npop)' %Muscle
'dC_adipose/dt = (1/vol_adipose).*(Q_adipose.*(C_arterial - C_adipose./Kpb_adipose)); C_adipose(0) = 0*ones(1,Npop)' %Adipose
'dC_bone/dt = (1/vol_bone).*(Q_bone.*(C_arterial - C_bone./Kpb_bone)); C_bone(0) = 0*ones(1,Npop)' %Bone
'dC_rest/dt = (1/vol_rest).*(Q_rest.*(C_arterial - C_rest./Kpb_rest)); C_rest(0) = 0*ones(1,Npop)' %Rest of body
'dC_venous/dt = (1/vol_venous).*(drug(t) + (Q_spleen+Q_gut+Q_liver).*C_liver./Kpb_liver + Q_kidney.*(C_vcd) + Q_brain.*(C_brain./Kpb_brain) + Q_heart.*(C_heart./Kpb_heart) + Q_skin.*(C_skin./Kpb_skin) + Q_muscle.*(C_muscle./Kpb_muscle) + Q_adipose.*(C_adipose./Kpb_adipose) + Q_bone.*(C_bone./Kpb_bone) + Q_rest.*(C_rest./Kpb_rest) - Q_venous.*C_venous); C_venous(0) = 0*ones(1,Npop)'%Venous compartment (collects all but gut and spleen lung)
%Liver (gets input from artery, gut, and spleen)
'dC_liver/dt = (1/vol_liver).*(Q_spleen.*(C_spleen./Kpb_spleen) + Q_gut.*(C_gut./Kpb_gut) + Q_liver.*C_arterial - (Q_spleen+Q_gut+Q_liver).*(C_liver./Kpb_liver) - clint_hep.*BW.*fub.*(C_liver./Kpb_liver)); C_liver(0) = 0*ones(1,Npop)'
%A separate compartment holds liver-metabolized drug
'dElim_liv/dt = (1./vol_liver).*clint_hep.*BW.*fub.*(C_liver./Kpb_liver)'
%Kidney
%This is the no-elimination model originally used in Tsamandouras et al., 2017
%'dC_kidney/dt = (Q_kidney.*(C_arterial - C_kidney./Kpb_kidney))./(Vol_kidney); C_kidney(0) = 0*ones(1,Npop)' %Temp kidney for now
%This is the four-part kidney model adapted from Scotcher with vascular and tubular sides
'dC_tpt/dt = (1./vol_t_pt).*(GFR.*fub.*C_arterial - Q_t_pt.*C_tpt - PkidAB.*papp_SF.*sa_t_pt.*(C_tpt)); C_tpt(0) = 0*ones(1,Npop)' %-MM(Ctpt) +
'dC_vpt/dt = (1./vol_v_pt).*(Q_kidney.*((1 - fub.*GFR./Q_kidney).*C_arterial - C_vpt) + PkidAB.*papp_SF.*sa_t_pt.*(C_tpt)); C_vpt(0) = 0*ones(1,Npop)' %MM(Ctpt) + %
'dC_tloh/dt = (1./vol_t_loh).*(Q_t_pt.*C_tpt - Q_t_loh.*C_tloh - PkidAB.*papp_SF.*sa_t_loh.*(C_tloh)); C_tloh(0) = 0*ones(1,Npop)'
'dC_vloh/dt = (1./vol_v_loh).*(Q_kidney.*(C_vpt - C_vloh) + PkidAB.*papp_SF.*sa_t_loh.*(C_tloh)); C_vloh(0) = 0*ones(1,Npop)' %
'dC_tdt/dt = (1./vol_t_dt).*(Q_t_loh.*C_tloh - Q_t_dt.*C_tdt - PkidAB.*papp_SF.*sa_t_dt.*(C_tdt)); C_tpt(0) = 0*ones(1,Npop)'
'dC_vdt/dt = (1./vol_v_dt).*(Q_kidney.*(C_vloh - C_vdt) + PkidAB.*papp_SF.*sa_t_dt.*(C_tdt)); C_vpt(0) = 0*ones(1,Npop)' %
'dC_tcd/dt = (1./vol_t_cd).*(Q_t_dt*C_tdt - Q_t_cd.*C_tcd - PkidAB.*papp_SF.*sa_t_cd.*(C_tcd)); C_tpt(0) = 0*ones(1,Npop)'
'dC_vcd/dt = (1./vol_v_cd).*(Q_kidney.*(C_vdt - C_vcd) + PkidAB.*papp_SF.*sa_t_cd.*(C_tcd)); C_vpt(0) = 0*ones(1,Npop)' %
%A separate compartment holds renally excreted drug
'dElim_kid/dt = (1./vol_t_cd).*(Q_t_cd.*C_tcd)'};
% This holds elements that are initialized inside the generated ODE file
% but whose values are replaced by the
P = {'PkidAB', PkidAB, 'PkidBA', PkidBA, 'Ht', Ht, 'fup', fup, 'logP',logP,...
'ionization',ionization,'pKa',pKa, 'dose', dose, 'Rbp', Rbp, 'clint_hep',...
clint_hep, 'infusion_duration', infusion_duration};
% DynaSim specification
spec = [];
spec.nodes.size = num_samples;
spec.nodes.equations = eqns;
spec.nodes.parameters = P;
% Run simulations
data = dsSimulate(spec,'tspan',tspan,'dt',.001,'save_data_flag',0,'solver','euler','precision','double');
%Plot profiles
%dsPlot(data);
t = data.time;
%Extract concentrations (ng/ml)
Cc = data.model_C_venous;
Cliv = data.model_C_liver;
Ccinit = data.model_C_venous(1);
%Extract kidney compartment concentrations
%Vascular kidney
Cvpt = data.model_C_vpt;
Cvloh = data.model_C_vloh;
Cvdt = data.model_C_vdt;
Cvcd = data.model_C_vcd;
%Tubular kidney
Ctpt = data.model_C_tpt;
Ctloh = data.model_C_tloh;
Ctdt = data.model_C_tdt;
Ctcd = data.model_C_tcd;
%Extract eliminated amounts
Elim_kid = data.model_Elim_kid;
Elim_liv = data.model_Elim_liv;
%Extract the remaining compartment concentrations (ng/ml)
C_venous = data.model_C_venous;
C_lung = data.model_C_lung;
C_arterial = data.model_C_arterial;
C_spleen = data.model_C_spleen;
C_gut = data.model_C_gut;
C_liver = data.model_C_liver;
C_brain = data.model_C_brain;
C_heart = data.model_C_heart;
C_skin = data.model_C_skin;
C_muscle = data.model_C_muscle;
C_adipose = data.model_C_adipose;
C_bone = data.model_C_bone;
C_rest = data.model_C_rest;
%C_kidney = data.model_C_kidney; %This is only for the transit-only kidney
%Obtain total drug mass for each compartment (ng)
M_venous = C_venous.*(data.model.fixed_variables.model_vol_venous);
M_lung = C_lung.*(data.model.fixed_variables.model_vol_lung);
M_arterial = C_arterial.*(data.model.fixed_variables.model_vol_arterial);
M_spleen = C_spleen.*(data.model.fixed_variables.model_vol_spleen);
M_gut = C_gut.*(data.model.fixed_variables.model_vol_gut);
M_liver = C_liver.*(data.model.fixed_variables.model_vol_liver);
M_brain = C_brain.*(data.model.fixed_variables.model_vol_brain);
M_heart = C_heart.*(data.model.fixed_variables.model_vol_heart);
M_skin = C_skin.*(data.model.fixed_variables.model_vol_skin);
M_muscle = C_muscle.*(data.model.fixed_variables.model_vol_muscle);
M_adipose = C_adipose.*(data.model.fixed_variables.model_vol_adipose);
M_bone = C_bone.*(data.model.fixed_variables.model_vol_bone);
M_rest = C_rest.*(data.model.fixed_variables.model_vol_rest);
M_vpt = Cvpt.*(data.model.fixed_variables.model_vol_t_pt);
M_vloh = Cvloh.*(data.model.fixed_variables.model_vol_t_loh);
M_vdt = Cvdt.*(data.model.fixed_variables.model_vol_t_dt);
M_vcd = Cvcd.*(data.model.fixed_variables.model_vol_t_cd);
%Tubular kidney
M_tpt = Ctpt.*(data.model.fixed_variables.model_vol_t_pt);
M_tloh = Ctloh.*(data.model.fixed_variables.model_vol_t_loh);
M_tdt = Ctdt.*(data.model.fixed_variables.model_vol_t_dt);
M_tcd = Ctcd.*(data.model.fixed_variables.model_vol_t_cd);
%Excreted mass (ng)
M_elim_kid = Elim_kid.*(data.model.fixed_variables.model_vol_t_cd);
M_elim_liv = Elim_liv.*(data.model.fixed_variables.model_vol_liver);
%M_kidney = C_kidney.*(data.model.fixed_variables.model_Vol_kidney);
%Create a total drug mass at all times from the above
Totdrug = (M_venous + M_lung + M_arterial ...
+ M_vpt + M_vloh + M_vdt + M_vcd ...
+ M_tpt + M_tloh + M_tdt + M_tcd ...
+ M_spleen + M_gut + M_liver ...
+ M_brain + M_heart + M_skin ...
+ M_muscle + M_adipose + M_bone + M_rest ...
+ M_elim_kid + M_elim_liv);
model_nca = noncompartmental_analysis(dose,C_venous,t);
if clinical_comparison_flag
data_nca = noncompartmental_analysis(dose,concentrations,sampling_times);
end
model_clr = M_elim_kid(end,:)./model_nca.AUC_obs;
model_clh = M_elim_liv(end,:)./model_nca.AUC_obs;
f_ren = model_clr./(model_clh+model_clr);
figure
subplot(2,1,1); semilogy(t./60,Cc); xlabel('time (h)'); ylabel('Log Venous Conc (ng/ml)'); %(xlim(tspan./60); ylim([10,max(Cc)*1.2]);)%
if clinical_comparison_flag
hold;
semilogy(sampling_times./60,concentrations, 'o'); xlim(tspan./60);
hold off;
end
subplot(2,1,2); semilogy(t./60,Cliv); xlabel('time (h)'); ylabel('Log Liver Conc (ng/ml)');
figure
subplot(2,2,1); plot(t./60,Cvpt); xlabel('time (h)'); ylabel('Vasc Proximal Tubule (ng/ml)');
subplot(2,2,2); plot(t./60,Cvloh); xlabel('time (h)'); ylabel('Vasc LoH (ng/ml)');
subplot(2,2,3); plot(t./60,Cvdt); xlabel('time (h)'); ylabel('Vasc Dist Tub (ng/ml)');
subplot(2,2,4); plot(t./60,Cvcd); xlabel('time (h)'); ylabel('Vasc Coll Duct (ng/ml)');
figure
subplot(2,2,1); plot(t./60,Ctpt); xlabel('time (h)'); ylabel('Proximal Tubule (ng/ml)');
subplot(2,2,2); plot(t./60,Ctloh); xlabel('time (h)'); ylabel('Loop of Henle (ng/ml)');
subplot(2,2,3); plot(t./60,Ctdt); xlabel('time (h)'); ylabel('Distal Tubule (ng/ml)');
subplot(2,2,4); plot(t./60,Ctcd); xlabel('time (h)'); ylabel('Collecting Duct (ng/ml)');
figure
plot(t,Totdrug); xlabel('time (h)'); ylabel('Total drug in system');
%% Analyze results
if num_samples==1, return; end
% Compute output metrics (frenal, CLr)
CLr = model_clr; %CLren; % [1 x nsims], TODO: get equations out of eqn object
frenal = f_ren; %PEU; % [1 x nsims], TODO: get equations out of eqn object
Cmax = max(C_arterial);
AUC = trapz(t,C_arterial);
model_name = data.model.specification.populations.name;
varied_params = cellfun(@(x)strrep(x,[model_name '_'],''),varied,'UniformOutput',false);
% Extract the varied parameters
P1 = data.model.parameters.(varied{1});
P = [P1'];
plab1=strrep(varied_params{1},'_','\_'); % [1 x nsims]
if length(varied)>1
P2 = data.model.parameters.(varied{2});
plab2=strrep(varied_params{2},'_','\_'); % [1 x nsims]
P = [P P2'];
Plabels = {plab1, plab2};
end
if length(varied)>2
P3 = data.model.parameters.(varied{3});
plab3=strrep(varied_params{3},'_','\_'); % [1 x nsims]
P = [P P3'];
Plabels = [Plabels, plab3];
end
if length(varied)>3
P4 = data.model.parameters.(varied{4});
plab4=strrep(varied_params{4},'_','\_'); % [1 x nsims]
P = [P P4'];
Plabels = [Plabels, plab4];
end
if length(varied)>4
P5 = data.model.parameters.(varied{5});
plab5=strrep(varied_params{5},'_','\_'); % [1 x nsims]
P = [P P5'];
Plabels = [Plabels, plab5];
end
% Sensitivity analysis
rho = partialcorr([CLr' frenal' Cmax' AUC'], P, zeros(num_samples,1));
Mlabels = {'CLr','frenal','Cmax', 'AUC'};
figure('position',[802 410 376 189]);
bar(rho'); ylabel('partial correlation')
set(gca,'xticklabel',Plabels);
legend(Mlabels{:},'Location','NorthEastOutside','Fontsize',12);
ylim([-1 1])
% Plot output metrics over parameter space
x = double(P1);
y = double(P2);
xgrid = linspace(min(x), max(x), 200);
ygrid = linspace(min(y), max(y), 200);
[X Y] = meshgrid(xgrid, ygrid);
z = double(CLr);
zlab = sprintf('Renal clearance');
Z = griddata(x,y,z, X,Y);
figure
surf(X,Y,Z,'EdgeAlpha',0); xlabel(plab1); ylabel(ylab); zlabel(zlab); title(zlab);
