Using Pumas-QSP: Global Sensitivity Analysis example

The package also enables users to consider additional model exploration based on Global Sensitivity Analysis (GSA).

Methodological background

This approach is helpful to further investigate how the model output reacts to variance in the model input. For example, we can analyse quantitatively how the variance in single input variables and how combinations of variances in several input variables affect output variables as shown in the figure below.

Influence of first input variable,
Influence of second input variable,
Influence of third input variable,
Influence of second and third input variable.

GSA objective

Here we present an example of performing GSA on the Hires model. Several global sensitivity analysis methods are supported, as listed here. For this example, we focus on the Sobol method.

Julia environment

First we prepare the environment by listening the packages we are using within the example.

# Packages used for this example
using PumasQSP
using Plots
using OrdinaryDiffEq
using ModelingToolkit
using DataFrames
using StatsPlots
using GlobalSensitivity
using Statistics

Model setup

The Hires model consists of eight states. We define the model here as an ODESystem.

# Defining QSP model
iv = @variables t # Independent variables
states = @variables y1(t) y2(t) y3(t) y4(t) y5(t) y6(t) y7(t) y8(t) # Dependent variables
ps = @parameters a b  # Parameters of the differential equations
D = Differential(t)
# The differential equations for the dependent variables
eqs = [ D(y1)  ~  (-a*y1 + 0.43*y2 + 8.32*y3 + 0.0007),
        D(y2)  ~  (a*y1 - 8.75*y2),
        D(y3)  ~  (-10.03*y3 + 0.43*y4 + 0.035*y5),
        D(y4)  ~  (8.32*y2 + a*y3 - 1.12*y4),
        D(y5)  ~  (-1.745*y5 + 0.43*y6 + 0.43*y7),
        D(y6)  ~  (-b*y6*y8 + 0.69*y4 + a*y5 - 0.43*y6 + 0.69*y7),
        D(y7)  ~  (b*y6*y8 - 1.81*y7),
        D(y8)  ~  (-b*y6*y8 + 1.81*y7)
]
# These are the values for the states and parameters
defs = Dict(y1 => 1.,
        y2 => 0.,
        y3 => 0.,
        y4 => 0.,
        y5 => 0.,
        y6 => 0.,
        y7 => 0.,
        y8 => 0.0057,
        a => 1.71,
        b => 280.0
)
# Model definition
@named model = ODESystem(eqs, t, states, ps; defaults = defs, checks = false)

Analysis

We solve the system for two fixed modeling setups in order to simulate data from two clinical conditions. We create again two Trial objects and implement a reduction where the solution of the ODE system is summarised for each state. Here this is the mean for the each state and they are calculated using the Statistics package.

# Clinical condition 1
trial1_tspan = (0.0, 10.)

# Data for trial 1
trail1_n = 10
trial1_saveat = range(trial1_tspan[1], trial1_tspan[2], length = trail1_n)
trial1_prob = ODEProblem(model, model.defaults, trial1_tspan, saveat = trial1_saveat)
trial1_sol = solve(trial1_prob, Tsit5(), reltol=1e-8, abstol=1e-8)
trial1_data = DataFrame(trial1_sol)

# Creating the Trial object for trial 1
trial1 = Trial(trial1_data, model;
            tspan = trial1_tspan,
            alg = Tsit5(),
            reduction = sol -> mean(sol, dims = 2)
	)

# Clinical condition 2
trial2_ps = [a => 1.61,
                b => 350.0]
trial2_tspan = (0.0, 20)

# Simulated data for this condition
trail2_n = 20
trial2_saveat = range(trial2_tspan[1], trial2_tspan[2], length = trail2_n)
trial2_prob = ODEProblem(model, model.defaults, trial2_tspan, trial2_ps, saveat = trial2_saveat)
trial2_sol = solve(trial2_prob, Tsit5(), reltol=1e-8, abstol=1e-8)
trial2_data = DataFrame(trial2_sol)

# Creating the Trial object for trial 2
trial2 = Trial(trial2_data, model;
                params = trial2_ps,
                tspan = trial2_tspan,
                alg = Tsit5(),
                reduction = sol -> mean(sol, dims = 2)
)
trials = [trial1, trial2]

Subsequently, we define the ranges for the parameters of the model to generate quasi-random samples using the Sobol sequence. Lastly, we define the QSPSensitivity object with the ranges of the parameters to be analyzed, and perform the gsa call.

# Create QSPSensitivity
sens = QSPSensitivity(model, trials, parameter_space = [
        a => (1.5, 1.8),
        b => (250, 300)
])
# Retrieve sensitivities
gsa_res = gsa(sens, Sobol(), batch = false, N = 1000)

Visualization

As we used the Sobol method with default settings, we can visualize the first and the total indices.

# Visualise GSA results
vals_first_a = gsa_res[1].S1[:,1]
vals_first_b = gsa_res[1].S1[:,2]
vals_total_a = gsa_res[1].ST[:,1]
vals_total_b = gsa_res[1].ST[:,2]

# Visualise GSA results
# First indices
p_first_a = bar(vals_first_a, grid = "off", ylabel = "Sensitivity (first)", xlabel = "States", label = "", title = "Parameter a")
p_total_a = bar(vals_total_a, grid = "off", ylabel = "Sensitivity (total)", xlabel = "States", label = "", title = "Parameter a")
plot(p_first_a, p_total_a, layout = (2,1))
# Total indices
p_first_b = bar(vals_first_b, grid = "off", ylabel = "Sensitivity (first)", xlabel = "States", label = "", title = "Parameter b")
p_total_b = bar(vals_total_b, grid = "off", ylabel = "Sensitivity (total)", xlabel = "States", label = "", title = "Parameter b")
plot(p_first_b, p_total_b, layout = (2,1))

Note: Instead of ranges, we can describe the respective area of the parameter space via design matrices, too. For details see GlobalSensitivity.jl.