Parameter Sensitivity

This notebook uses the direct sensitivity solver in AutoReduce for a one-state decay model. Local sensitivity analysis is performed by solving the sensitivity equations alongside the original system, then comparing the result against an analytical expression.

Model

The model is

\[\dot{x} = -kx,\]

with solution \(x(t)=x_0e^{-kt}\). The sensitivity of the state with respect to \(k\) is

\[\frac{\partial x}{\partial k} = -t x_0 e^{-kt}.\]

[ ]:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sympy import Symbol

from autoreduce import System, solve_sensitivity

x = Symbol("x")
k = Symbol("k")

system = System(
    [x],
    [-k * x],
    params=[k],
    params_values=[0.4],
    x_init=[5.0],
    C=np.array([[1.0]]),
)

Compute local sensitivities

solve_sensitivity returns an array indexed by time point, parameter, and state.

[ ]:

timepoints = np.linspace(0.0, 5.0, 1024)
sensitivities = solve_sensitivity(system, timepoints, normalize=False)

sensitivities.shape

Compare against analytical coefficients

For this model, the analytical expression is available, so it provides a compact check on the numerical sensitivity calculation.

[ ]:

x0 = system.x_init[0]
k_value = system.params_values[0]
analytical = -timepoints * x0 * np.exp(-k_value * timepoints)
numerical = sensitivities[:, 0, 0]

comparison = np.column_stack([timepoints, numerical, analytical])
comparison[:8]

Verify the sensitivity calculation

The maximum absolute error gives a direct numerical check of the computed sensitivity against the analytical expression.

[ ]:

absolute_error = np.abs(numerical - analytical)
max_absolute_error = float(np.max(absolute_error))
tolerance = 1e-6

assert max_absolute_error < tolerance
{
    "max_absolute_error": max_absolute_error,
    "tolerance": tolerance,
    "passes": max_absolute_error < tolerance,
}

Plot a sensitivity heatmap

The heatmap shows how the local sensitivity of \(x\) to \(k\) changes over time.

[ ]:

fig, ax = plt.subplots(figsize=(10, 2.2))
sns.heatmap(
    numerical.reshape(1, -1),
    cmap="vlag",
    center=0,
    xticklabels=128,
    yticklabels=[r"$\partial x / \partial k$"],
    cbar_kws={"label": "Sensitivity"},
    ax=ax,
)
ax.set_xlabel("Time index")
ax.set_ylabel("State/parameter pair")
ax.set_title("Sensitivity of x to k over time")
plt.tight_layout()