# 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,

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,

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.