Global Sensitivity Analysis on HCV model
2021-04-26
using Pumas, GlobalSensitivity, CairoMakie, PumasPlots
Introduction
In this tutorial, we will cover running global sensitivity analysis on the HCV model.
Model Code
peg_inf_model = @model begin @param begin logKa ∈ RealDomain() logKe ∈ RealDomain() logVd ∈ RealDomain() logn ∈ RealDomain() logd ∈ RealDomain() logc ∈ RealDomain() logEC50 ∈ RealDomain() ω²Ka ∈ RealDomain(lower=0.0) ω²Ke ∈ RealDomain(lower=0.0) ω²Vd ∈ RealDomain(lower=0.0) ω²n ∈ RealDomain(lower=0.0) ω²δ ∈ RealDomain(lower=0.0) ω²c ∈ RealDomain(lower=0.0) ω²EC50 ∈ RealDomain(lower=0.0) σ²PK ∈ RealDomain(lower=0.0) σ²PD ∈ RealDomain(lower=0.0) end @random begin ηKa ~ Normal(0.0, sqrt(ω²Ka)) ηKe ~ Normal(0.0, sqrt(ω²Ke)) ηVd ~ Normal(0.0, sqrt(ω²Vd)) ηn ~ Normal(0.0, sqrt(ω²n)) ηδ ~ Normal(0.0, sqrt(ω²δ)) ηc ~ Normal(0.0, sqrt(ω²c)) ηEC50 ~ Normal(0.0, sqrt(ω²EC50)) end @pre begin # constants p = 100.0 d = 0.001 e = 1e-7 s = 20000.0 log_Ka = logKa + ηKa log_Ke = logKe + ηKe log_Vd = logVd + ηVd log_n = logn + ηn log_d = logd + ηδ log_c = logc + ηc log_EC50 = logEC50 + ηEC50 end @init begin T = exp(log_c + log_d)/(p*e) I = (s*e*p - d*exp(log_c + log_d))/(p*exp(log_d)*e) W = (s*e*p - d*exp(log_c + log_d))/(exp(log_c + log_d)*e) end @dynamics begin X' = -exp(log_Ka)*X A' = exp(log_Ka)*X - exp(log_Ke)*A T' = s - T*(e*W + d) I' = e*W*T - exp(log_d)*I W' = p/((A/exp(log_Vd)/exp(log_EC50) + 1e-100)^exp(log_n) + 1)*I - exp(log_c)*W end @derived begin conc = @. A/exp(log_Vd) log10W = @. log10(max.(W, 1e-5)) yPK ~ @. Normal(A/exp(log_Vd), sqrt(σ²PK)) yPD ~ @. Normal(log10W, sqrt(σ²PD)) nca1 := @nca conc cmax_pk = NCA.cmax(nca1) nca2 := @nca log10W cmax_pd = NCA.cmax(nca2) end end
PumasModel Parameters: logKa, logKe, logVd, logn, logd, logc, logEC50, ω²Ka, ω²Ke, ω ²Vd, ω²n, ω²δ, ω²c, ω²EC50, σ²PK, σ²PD Random effects: ηKa, ηKe, ηVd, ηn, ηδ, ηc, ηEC50 Covariates: Dynamical variables: X, A, T, I, W Derived: conc, log10W, yPK, yPD, cmax_pk, cmax_pd Observed: conc, log10W, yPK, yPD, cmax_pk, cmax_pd
We generate a population to simulate with.
peg_inf_dr = DosageRegimen(180.0) t = collect(0.0:1.0:28.0) _pop = map(i -> Subject(id=i, observations=(yPK=[], yPD=[]), events=peg_inf_dr, time=t), 1:3)
Population Subjects: 3 Observations: yPK, yPD
GSA
Let's define the parameters now, we'll fix the ω parameters to neutralize effect of random effects and run the GSA only on the population parameters logKa ,logKe ,logVd ,logn ,logd ,logc and logEC50.
param_PKPD = ( logKa = log(0.80), logKe = log(0.15), logVd = log(100.0), logn = log(2.0), logd = log(0.20), logc = log(7.0), logEC50 = log(0.12), # random effects variance parameters, set to 0 to neutralize effect of random effects ω²Ka = 0.0, ω²Ke = 0.0, ω²Vd = 0.0, ω²n = 0.0, ω²δ = 0.0, ω²c = 0.0, ω²EC50 = 0.0, # variance parameter in proportional error model σ²PK = 0.04, σ²PD = 0.04)
(logKa = -0.2231435513142097, logKe = -1.8971199848858813, logVd = 4.605170 185988092, logn = 0.6931471805599453, logd = -1.6094379124341003, logc = 1. 9459101490553132, logEC50 = -2.120263536200091, ω²Ka = 0.0, ω²Ke = 0.0, ω²V d = 0.0, ω²n = 0.0, ω²δ = 0.0, ω²c = 0.0, ω²EC50 = 0.0, σ²PK = 0.04, σ²PD = 0.04)
We'll run the Morris method and Sobol method. The derived variables cmax_pk and cmax_pd which are CMAX outputs of NCA will be analysed.
Let's take a look at the simulation of the model to ensure everything is working as expected. We'll use threading across subjects for simulation and also utlize this parallelism down the line in GSA as well.
simdata = simobs(peg_inf_model, _pop, param_PKPD, ensemblealg=EnsembleThreads()) sim_plot(peg_inf_model ,simdata, observations=[:yPD,:yPK], all = true)
To run the GSA we'll define the parameter ranges for our parameters of interest. Let's set the range between half and twice of the initial values.
p_range_low = (logKa = log(0.80)/2, logKe = log(0.15)/2, logVd = log(100.0)/2, logn = log(2.0)/2, logd = log(0.20)/2, logc = log(7.0)/2, logEC50 = log(0.12)/2, ) p_range_high = (logKa = log(0.80)*2, logKe = log(0.15)*2, logVd = log(100.0)*2, logn = log(2.0)*2, logd = log(0.20)*2, logc = log(7.0)*2, logEC50 = log(0.12)*2, )
(logKa = -0.4462871026284194, logKe = -3.7942399697717626, logVd = 9.210340 371976184, logn = 1.3862943611198906, logd = -3.2188758248682006, logc = 3. 8918202981106265, logEC50 = -4.240527072400182)
Now, we are ready to run GSA on our model.
The Morris Method.
morris_ = Pumas.gsa( peg_inf_model, _pop,param_PKPD, GlobalSensitivity.Morris(), [:cmax_pk,:cmax_pd], p_range_low, p_range_high, ensemblealg=EnsembleThreads())
Morris Sensitivity Analysis
Means (μ)
2×8 DataFrame
Row │ dv_name logKa logKe logVd logn logd
⋯
│ Any Float64 Float64 Float64 Float64 Float64
⋯
─────┼─────────────────────────────────────────────────────────────────────
─────
1 │ cmax_pk 0.215144 -0.658552 -2.16182 1.43616e-8 -3.46422
e-9 ⋯
2 │ cmax_pd 0.000326863 0.00135321 -0.0205303 -0.011866 -0.42875
2
2 columns om
itted
Means star (μ*)
2×8 DataFrame
Row │ dv_name logKa logKe logVd logn logd
lo ⋯
│ Any Float64 Float64 Float64 Float64 Float64
Fl ⋯
─────┼─────────────────────────────────────────────────────────────────────
─────
1 │ cmax_pk 0.215144 0.658552 2.16182 5.42518e-8 3.86123e-9
3. ⋯
2 │ cmax_pd 0.000326863 0.00181066 0.0205303 0.0156105 0.428752
0.
2 columns om
itted
Variances
2×8 DataFrame
Row │ dv_name logKa logKe logVd logn logd
l ⋯
│ Any Float64 Float64 Float64 Float64 Float64
F ⋯
─────┼─────────────────────────────────────────────────────────────────────
─────
1 │ cmax_pk 0.0591453 0.416196 6.01715 8.91614e-15 3.50687e-17
1 ⋯
2 │ cmax_pd 1.42465e-7 1.8728e-5 0.00131828 0.00110238 0.00148458
0
2 columns om
itted
Let's create bar plots of the means of Morris Method output to visualize the output.
keys_ = [keys(p_range_low)...] cmax_pk_meansstar = [morris_.means_star[1,:][key] for key in keys_] cmax_pd_meansstar = [morris_.means_star[2,:][key] for key in keys_] fig = Figure(resolution = (1200, 400)) plot_pk = barplot(fig[1,1], string.(keys_), cmax_pk_meansstar, axis = (xlabel ="Parameters", title="Cmax PK Morris means")) plot_pd = barplot(fig[1,2], string.(keys_), cmax_pd_meansstar, axis = (xlabel ="Parameters", title="Cmax PD Morris means")) display(fig)
We observe that PK Cmax is most sensitive to logVd, logKa and logKe parameters where as PD Cmax is most sensitive to logc, logd and logKe.
The Sobol Method
sobol_ = Pumas.gsa( peg_inf_model, _pop, param_PKPD, GlobalSensitivity.Sobol(), [:cmax_pk,:cmax_pd], p_range_low, p_range_high, N=4000, ensemblealg=EnsembleThreads())
Sobol Sensitivity Analysis
First Order Indices
2×8 DataFrame
Row │ dv_name logKa logKe logVd logn logd
l ⋯
│ Any Float64 Float64 Float64 Float64 Float64
F ⋯
─────┼─────────────────────────────────────────────────────────────────────
─────
1 │ cmax_pk 0.00019134 0.0103614 0.961884 7.19654e-9 -1.33176e-7
- ⋯
2 │ cmax_pd -0.00111712 -0.095498 0.0176666 0.00241971 0.368883
2 columns om
itted
Total Order Indices
2×8 DataFrame
Row │ dv_name logKa logKe logVd logn logd
log ⋯
│ Any Float64 Float64 Float64 Float64 Float64
Flo ⋯
─────┼─────────────────────────────────────────────────────────────────────
─────
1 │ cmax_pk 0.000469732 0.037356 0.988439 -1.34007e-7 1.98202e-6
1.0 ⋯
2 │ cmax_pd 0.000810145 0.127653 0.0121693 -0.00256567 0.318561
0.6
2 columns om
itted
We use heatmaps to visualize the Sobol output
cmax_pk_s1 = [sobol_.first_order[1,:][key] for key in keys_] cmax_pd_s1 = [sobol_.first_order[2,:][key] for key in keys_] cmax_s1 = hcat(cmax_pk_s1, cmax_pd_s1) fig = Figure(resolution = (1200, 400)) plot_s1 = heatmap(fig[1,1], 1:2, 1:7, cmax_s1, axis = (xticks = (1:2, ["Cmax_PK","Cmax_PD"]), yticks = (1:7, string.(keys_)), title="First Order Indices"), colormap = "darkrainbow") cmax_pk_st = [sobol_.total_order[1,:][key] for key in keys_] cmax_pd_st = [sobol_.total_order[2,:][key] for key in keys_] cmax_st = hcat(cmax_pk_st, cmax_pd_st) plot_st = heatmap(fig[1,2], 1:2, 1:7, cmax_st, axis = (xticks = (1:2, ["Cmax_PK","Cmax_PD"]), yticks = (1:7, string.(keys_)), title="First Order Indices"), colormap = "darkrainbow") display(fig)
Conclusion
For PK Cmax logVd comes out to be the most dominant parameter by a large margin as evident by both First and Total Order indices of Sobol. PD Cmax is most sensitive to logc and logd and the Total Order indices show logKe is also important, hence we can deduce that logKe has a greater influence due to its interactions as compared to its contribution individually.