Clarification regarding Sensitivity Analysis of a coupled ODE system using the odesensitivity function

I need to confirm whether the the way I wrote the code for the sensitivity analysis of an ODE system is correct or not. I have done the sensitivity analysis using the odeSensitivity introduced in R2024a of MATLAB. I also want to know what the values exactly mean regarding the sensitivity of the system. The equations for the system are as follows

clear all; clc;

% Define Parameters
params = [1,                 % S_in (mgC/L)
          1 / (1.67e-5),     % tau (Residence time in sec)
          1e-4,              % mu_max (Max bacterial growth rate in s^-1)
          0.1,               % K_S (Half-saturation constant for substrate in mgC/L)
          0.6,               % Y_B (Bacterial yield coefficient)
          1e-5,              % lambda_B (Decay rate constant for bacteria in s^-1)
          1e-3,              % mu_max_graz (Max grazer growth rate in s^-1)
          10,                % K_bac (Half-saturation constant for grazers in mgC/L)
          1.0,               % Y_G (Grazer yield coefficient)
          5e-6];             % lambda_G (Decay rate constant for grazers in s^-1)

% Initial Conditions
S0 = 0; % Substrate concentration (mgC/L)
B0 = 0.01;   % Bacteria concentration (mgC/L)
G0 = 0.01;   % Grazer concentration (mgC/L)
initial_conditions = [S0, B0, G0];

% Time Span
tspan = [0, 4.32e6]; % 50 days in seconds

% Define odeSensitivity Object
sensitivity_obj = odeSensitivity();

% Create the ode object
F = ode(ODEFcn=@(t, C, p) retentostat_odes(t, C, p), ...
        InitialValue=initial_conditions, ...
        Parameters=params, ...
        Sensitivity=sensitivity_obj);

% Solve the ODE system
S = solve(F, tspan(1), tspan(2)); % Solve for the time span (50 days)

% Extract the baseline solution
S_baseline = S.Solution; % State variables at the baseline solution (initial params)

% Extract the sensitivity results
S_sensitivities = S.Sensitivity; % Sensitivity results for each parameter

% Visualizing Sensitivities
param_names = {'c_{in}', '\tau', '\mu_{max}^{bac}', 'K_{sub}', 'Y_{bac}', ...
               '\mu_{dec}^{bac}', '\mu_{max}^{graz}', 'K_{bac}', ...
               'Y_{graz}', '\mu_{dec}^{graz}'};

% Plot Sensitivity Results
figure;
% Set figure size to 1920x1080 pixels
set(gcf, 'Units', 'pixels', 'Position', [100, 100, 1920, 1080]);
for i = 1:length(param_names)
    subplot(2, ceil(length(param_names) / 2), i);
    hold on;
    for j = 1:3 % For each state variable (S, B, G)
        plot(S.Time / (3600 * 24), squeeze(S_sensitivities(j, i, :)),'Linewidth',1.25); % Time in days
    end
    title(param_names{i});
    xlabel('Time (days)', 'FontSize', 12, 'FontWeight', 'bold');
    ylabel('Sensitivity', 'FontSize', 12, 'FontWeight', 'bold');
    
    hold off;
    ay2 = gca;  
    ay2.FontSize = 12; % Set font size of ticks
    ay2.FontWeight = 'bold'; % Set font weight of ticks to bold
end
sgtitle('Parameter Sensitivities for Chemostat System', 'FontSize', 12, 'FontWeight', 'bold');
h = legend({'c_{sub}', 'c_{bac}', 'c_{graz}'});
set(h, 'Position', [0.92, 0.47, 0.01, 0.01]);
drawnow;

% Save the figure with 4k resolution
print('chemo_senst.png', '-dpng', '-r600');

% ODE function
function dCdt = retentostat_odes(t, C, p)
    % Extract state variables
    S = C(1); % Substrate concentration
    B = C(2); % Bacteria concentration
    G = C(3); % Grazer concentration
    
    % Reaction rates using indexed parameters
    sub_rate = ((p(3) / p(5)) * (S / (p(4) + S)) * B);   % Substrate utilization rate
    grazing_rate = G * (p(7) / p(9)) * (B / (p(8) + B));  % Grazing rate
    
    % ODEs
    dSdt = (1 / p(2)) * (p(1) - S) - sub_rate;  % Substrate dynamics
    dBdt = p(5) * sub_rate - grazing_rate - p(6) * B - B/p(2);  % Bacteria dynamics
    dGdt = p(9) * grazing_rate - p(10) * G - G/p(2);  % Grazer dynamics

    % Return derivatives
    dCdt = [dSdt; dBdt; dGdt];
end

Answer 1

[Torsten]

Start with the simple example

y'(t) = p(1)*y(t), y(0) = 1.

The solution is

y(t) = exp(p(1)*t).

The (local) sensitivity is defined as

sensitivity(t) = dy/d(p(1)) 

i.e. it is an indicator of how "sensible" y will react if parameter p(1) is changed at time t ( = delta y / delta p).

For the example above, you get

dy/d(p(1)) = d/d(p(1)) (exp(p(1)*t)) = t*exp(p(1)*t)

- and both curves agree in the below graphics.

I suggest you turn the initial condition y0 = 1 into a second parameter and analyze this extended example.

% Define Parameters
params = 1;
% Initial Conditions
y0 = 1; 
initial_conditions = y0;
% Time Span
tspan = [0, 5]; 
% Define odeSensitivity Object
sensitivity_obj = odeSensitivity();
% Create the ode object
F = ode(ODEFcn=@(t, y, p) odes(t, y, p), ...
        InitialValue=initial_conditions, ...
        Parameters=params, ...
        Sensitivity=sensitivity_obj);
% Solve the ODE system
S = solve(F, tspan(1), tspan(2)); % Solve for the time span (50 days)
% Extract the baseline solution
S_baseline = S.Solution; % State variables at the baseline solution (initial params)
% Extract the sensitivity results
S_sensitivities = S.Sensitivity; % Sensitivity results for each parameter
% Visualizing Sensitivities
param_names = {'a'};
% Plot Sensitivity Results
figure
hold on
plot(S.Time , squeeze(S_sensitivities(1,1, :)))
plot(S.Time , S.Time.*exp(params*S.Time))
hold off
title(param_names{1});
xlabel('Time');
ylabel('Sensitivity');
grid on

●

% ODE function
function dydt = odes(t, y, p)
    dydt = p(1)*y;
end