Sensitivity Analysis

Sensitivity Analysis#

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

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.
Running the simulation for 100 sets of parameters sampled from the parameter space.
Using Autoemulate to find the best emulator for this simulation
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 :

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).
L (Inductance): Represents the inertial effects of blood flow, capturing how blood resists changes in its velocity (Analogous to electrical inductor).
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#

Neumann boundary condition : Specifies the derivative of the variable at the boundary.
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.

\(\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", "GP"], log_level="error")  # remove models argument to use all models
best = ae.best_result()
print(best.model_name)

GaussianProcess