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.000578 0.002691
1 pressure td S1 0.013941 0.011895
2 pressure amp S1 0.925220 0.064498
3 pressure dt S1 0.000088 0.001370
4 pressure C S1 0.009650 0.007560
5 pressure R S1 0.036995 0.018027
6 pressure L S1 0.001157 0.002949
7 pressure R_o S1 0.000281 0.001064
8 pressure p_o S1 0.000836 0.000689
0 pressure T ST 0.001020 0.000140
1 pressure td ST 0.019115 0.002527
2 pressure amp ST 0.936366 0.061310
3 pressure dt ST 0.000303 0.000043
4 pressure C ST 0.010006 0.001353
5 pressure R ST 0.043748 0.006191
6 pressure L ST 0.001332 0.000127
7 pressure R_o ST 0.000137 0.000016
8 pressure p_o ST 0.000055 0.000009
0 pressure (T, td) S2 0.000098 0.004150
1 pressure (T, amp) S2 0.000793 0.005033
2 pressure (T, dt) S2 0.000249 0.004215
3 pressure (T, C) S2 0.000186 0.004289
4 pressure (T, R) S2 0.000058 0.004261
5 pressure (T, L) S2 0.000279 0.004228
6 pressure (T, R_o) S2 0.000265 0.004220
7 pressure (T, p_o) S2 0.000264 0.004228
8 pressure (td, amp) S2 0.004137 0.019881
9 pressure (td, dt) S2 0.001094 0.017278
10 pressure (td, C) S2 0.001349 0.017142
11 pressure (td, R) S2 0.001443 0.017415
12 pressure (td, L) S2 0.001018 0.017186
13 pressure (td, R_o) S2 0.001128 0.017250
14 pressure (td, p_o) S2 0.001147 0.017285
15 pressure (amp, dt) S2 -0.000533 0.071528
16 pressure (amp, C) S2 -0.000726 0.071934
17 pressure (amp, R) S2 0.005201 0.080206
18 pressure (amp, L) S2 -0.000557 0.072485
19 pressure (amp, R_o) S2 -0.000923 0.071695
20 pressure (amp, p_o) S2 -0.001516 0.071795
21 pressure (dt, C) S2 -0.000141 0.002582
22 pressure (dt, R) S2 0.000034 0.002598
23 pressure (dt, L) S2 -0.000098 0.002579
24 pressure (dt, R_o) S2 -0.000141 0.002575
25 pressure (dt, p_o) S2 -0.000131 0.002574
26 pressure (C, R) S2 -0.001490 0.009786
27 pressure (C, L) S2 -0.000956 0.010500
28 pressure (C, R_o) S2 -0.000943 0.010492
29 pressure (C, p_o) S2 -0.000905 0.010521
30 pressure (R, L) S2 -0.001111 0.029082
31 pressure (R, R_o) S2 -0.001302 0.029004
32 pressure (R, p_o) S2 -0.001097 0.029011
33 pressure (L, R_o) S2 0.000082 0.004344
34 pressure (L, p_o) S2 0.000144 0.004344
35 pressure (R_o, p_o) S2 0.000236 0.001956
sa.plot_sobol(sobol_df, index="ST", figsize=figsize) 
../../_images/af15591124f6213c3abc1d4f0dcb87878421360765df401d05958e38ac53b43a.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 -16.771048 21.454834 21.596975 1.066423
1 pressure td 79.862778 83.085632 69.878227 3.799350
2 pressure amp 694.629211 694.629211 99.040108 5.858261
3 pressure dt -4.638096 13.009188 15.113632 0.526546
4 pressure C -63.811440 64.648178 40.131725 2.096304
5 pressure R -132.567734 133.197311 87.748726 5.516351
6 pressure L 24.545521 24.545521 10.844708 0.698015
7 pressure R_o -2.794015 8.773822 9.840091 0.310575
8 pressure p_o -3.208048 5.181366 5.456545 0.196669
sa.plot_morris(morris_df, figsize=figsize)
../../_images/acf6017e4d0258d810d65c0252a8729fe14b9c60f6b21a53688409a3e78c2813.png