Sensitivity Analysis#

In this tutorial we demonstrate how to perform sensitivity analysis as part of the AutoEmulate workflow. The tutorial covers:

  1. Setting up an example simulation: here we use our “FlowProblem” simulator. This is a cardiovascular modelling example, simulating a blood vessel divided into 10 compartments. This allows for the study of the pressure and flow rate at various points in the tube. See “The Flow Problem” below for more details.

  2. Running the simulation for 100 sets of parameters sampled from the parameter space.

  3. Using Autoemulate to find the best emulator for this simulation

  4. Performing sensitivity analysis.

The Flow Problem

In the field of cardiovascular modeling, capturing the dynamics of blood flow and the associated pressures and volumes within the vascular system is crucial for understanding heart function and disease. This simulator simulates a vessel divided to 10 compartments.

Parameters#

The simulation parameters include :

  1. R (Resistance): Represents the resistance to blood flow in blood vessels, akin to the hydraulic resistance caused by vessel diameter and blood viscosity (Analogous to electrical resistor).

  2. L (Inductance): Represents the inertial effects of blood flow, capturing how blood resists changes in its velocity (Analogous to electrical inductor).

  3. C (Capacitance): Represents the compliance or elasticity of blood vessels, primarily large arteries, which store and release blood volume with changes in pressure (Analogous to a capacitor).

Boundary conditions#

  1. Neumann boundary condition : Specifies the derivative of the variable at the boundary.

  2. Dirichlet Boundary condition : Specifies the value of the variable directly at the boundary.

The setup#

The input flow rate in each compartment is \(Q_i(t)\) for the \(i^{th}\) compartment and the output flow rate is \(Q_{i+1}(t)\).

\(Q_0(t) = \begin{cases} A \cdot \sin^2\left(\frac{\pi}{t_d} t\right), & \text{if } 0 \leq t < T, \\ 0, & \text{otherwise}. \end{cases}\)

Where:

  • \(Q_0(t)\) is the input pulse function (flow rate) at time t

  • A is the amplitude of the pulse

  • \(t_d\) is the pulse duration.

Solve#

Pressure in Each Compartment (\(P_i\)): This determines how the pressure in each compartment evolves over time, based on the inflow (\(Q_i(t)\)) and the outflow (\(Q_{i+1}(t)\)). where \(i\) is the number of compartment.

Circuit Diagram

\(\frac{dP_i}{dt} = \frac{1}{C_n} \left( Q_i(t) - Q_{i+1}(t) \right)\) where, \(C_n = \frac{C}{n_\text{comp}}\)

Flow rate equation (\(Q_i\)): This governs how the flow in each compartment changes over time, depending on the pressures in the neighboring compartments and the resistance and inertance properties of each compartment.

\(\frac{dQ_i}{dt} = \frac{1}{L_n} \left( P_i - P_{10} - R_n Q_i(t) \right)\), where \(L_n = \frac{L}{n_\text{comp}}, \quad R_n = \frac{R}{n_\text{comp}}\)

from autoemulate.core.compare import AutoEmulate
from autoemulate.core.sensitivity_analysis import SensitivityAnalysis
from autoemulate.simulations.flow_problem import FlowProblem
import warnings
warnings.filterwarnings("ignore")
figsize = (9, 5)

Set up the simulation parameters and ranges:

parameters_range = {
    "T": (0.5, 2.0), # Cardiac cycle period (s)
    "td": (0.1, 0.5), # Pulse duration (s)
    "amp": (100.0, 1000.0), # Amplitude (e.g., pressure or flow rate)
    "dt": (0.0001, 0.01), # Time step (s)
    "C": (20.0, 60.0), # Compliance (unit varies based on context)
    "R": (0.01, 0.1), # Resistance (unit varies based on context)
    "L": (0.001, 0.005), # Inductance (unit varies based on context)
    "R_o": (0.01, 0.05), # Outflow resistance (unit varies based on context)
    "p_o": (5.0, 15.0) # Initial pressure (unit varies based on context)
}
output_names = ["pressure"]

simulator = FlowProblem(
    parameters_range=parameters_range,
    output_names=output_names,
    log_level="error"
)

Run the simulation for 100 sets of parameters sampled from the parameter space:

x = simulator.sample_inputs(100)
y, _ = simulator.forward_batch(x)
print(x.shape, y.shape)
torch.Size([100, 9]) torch.Size([100, 1])

Use AutoEmulate to find the best emulator for this simulation:

ae = AutoEmulate(x, y, models=["MLP", "GaussianProcessRBF"], log_level="error")  # remove models argument to use all models
best = ae.best_result()
print(best.model_name)
GaussianProcessRBF

Sensitivity Analysis#

  1. Define the problem by creating a dictionary which contains the names and the boundaries of the parameters

  2. Evaluate the contribution of each parameter via the Sobol and Morris methods.

problem = {
    'num_vars': simulator.in_dim,
    'names': simulator.param_names,
    'bounds': simulator.param_bounds,
    'output_names': simulator.output_names,
}
sa = SensitivityAnalysis(best.model, problem=problem)

Sobol metrics:

  • \(S_1\): First-order sensitivity index.

  • \(S_2\): Second-order sensitivity index.

  • \(S_t\): Total sensitivity index.

Sobol interpretation:

  • \(S_1\) values sum to ≤ 1.0 (exact fraction of variance explained)

  • \(S_t - S_1\) = interaction effects involving that parameter

  • Large \(S_t - S_1\) gap indicates strong interactions

sobol_df = sa.run("sobol")
sobol_df
output parameter index value confidence
0 pressure T S1 0.000282 0.002487
1 pressure td S1 0.009502 0.010796
2 pressure amp S1 0.933604 0.075751
3 pressure dt S1 -0.000147 0.000949
4 pressure C S1 0.007805 0.008617
5 pressure R S1 0.031847 0.019118
6 pressure L S1 0.000718 0.004544
7 pressure R_o S1 -0.001959 0.001884
8 pressure p_o S1 0.001078 0.001064
0 pressure T ST 0.000737 0.000090
1 pressure td ST 0.013830 0.001611
2 pressure amp ST 0.946087 0.060683
3 pressure dt ST 0.000133 0.000020
4 pressure C ST 0.010015 0.001147
5 pressure R ST 0.037964 0.004165
6 pressure L ST 0.002349 0.000311
7 pressure R_o ST 0.000420 0.000059
8 pressure p_o ST 0.000134 0.000021
0 pressure (T, td) S2 -0.000196 0.003723
1 pressure (T, amp) S2 0.000296 0.004163
2 pressure (T, dt) S2 -0.000139 0.003687
3 pressure (T, C) S2 -0.000096 0.003707
4 pressure (T, R) S2 0.000104 0.003704
5 pressure (T, L) S2 -0.000131 0.003703
6 pressure (T, R_o) S2 -0.000074 0.003670
7 pressure (T, p_o) S2 -0.000136 0.003691
8 pressure (td, amp) S2 0.003465 0.020166
9 pressure (td, dt) S2 0.001047 0.015953
10 pressure (td, C) S2 0.001428 0.016015
11 pressure (td, R) S2 0.001261 0.015833
12 pressure (td, L) S2 0.001153 0.016092
13 pressure (td, R_o) S2 0.001430 0.015879
14 pressure (td, p_o) S2 0.001057 0.015961
15 pressure (amp, dt) S2 -0.000487 0.082832
16 pressure (amp, C) S2 0.000748 0.083603
17 pressure (amp, R) S2 0.005156 0.088109
18 pressure (amp, L) S2 0.000989 0.083425
19 pressure (amp, R_o) S2 0.001225 0.082047
20 pressure (amp, p_o) S2 -0.002019 0.082637
21 pressure (dt, C) S2 0.000176 0.001601
22 pressure (dt, R) S2 0.000110 0.001606
23 pressure (dt, L) S2 0.000145 0.001582
24 pressure (dt, R_o) S2 0.000156 0.001579
25 pressure (dt, p_o) S2 0.000125 0.001583
26 pressure (C, R) S2 0.000924 0.012064
27 pressure (C, L) S2 0.000937 0.011624
28 pressure (C, R_o) S2 0.001004 0.011551
29 pressure (C, p_o) S2 0.000949 0.011606
30 pressure (R, L) S2 -0.000160 0.026177
31 pressure (R, R_o) S2 -0.000016 0.025974
32 pressure (R, p_o) S2 0.000026 0.026026
33 pressure (L, R_o) S2 0.000041 0.007464
34 pressure (L, p_o) S2 -0.000053 0.007460
35 pressure (R_o, p_o) S2 0.001588 0.002513
sa.plot_sobol(sobol_df, index="ST", figsize=figsize) 
../../_images/cd0f1ca45551872f4400db97a21b3a77e3f2fbb033c835754cc9c062c959908d.png

You can also save the plot directly to a file by passing the fname argument to the plotting function.

sa.plot_sobol(sobol_df, index="ST", figsize=figsize, fname="sobol_plot.png") 

Morris Interpretation:

  • High \(\mu^*\), Low \(\sigma\): Important parameter with linear/monotonic effects

  • High \(\mu^*\), High \(\sigma\): Important parameter with non-linear effects or interactions

  • Low \(\mu^*\), High \(\sigma\): Parameter involved in interactions but not individually important

  • Low \(\mu^*\), Low \(\sigma\): Unimportant parameter

morris_df = sa.run("morris")
morris_df
output parameter mu mu_star sigma mu_star_conf
0 pressure T -13.390306 17.969063 17.221891 0.696981
1 pressure td 72.336555 74.324181 58.860565 3.594071
2 pressure amp 674.267944 674.267944 91.172890 5.842741
3 pressure dt 1.585763 8.564500 10.739902 0.417649
4 pressure C -60.633743 62.018288 41.782467 2.625212
5 pressure R -122.803902 123.726814 71.418617 4.167825
6 pressure L 20.448790 29.310217 30.362619 1.486954
7 pressure R_o -2.793113 15.379130 18.780693 0.681314
8 pressure p_o -1.855209 8.475296 10.467030 0.365496
sa.plot_morris(morris_df, figsize=figsize)
../../_images/d93878dafff1087ecf792ef37dc34d67f09365584fbeecd7ca85d49a9ddc2e2d.png