SENSITIVITY ANALYSIS OF A SYSTEM OF AN ODE USING NORMALIZED SENSITIVITY FUNCTION

I want to perform the sensitivity analysis on the parameters of an ODE SIR model using normalized sensitivity function. I used the following code below. When it was run, I got this error: "undefined functionnor variable 'p', error in epidemic1 line8, F=ode45(@epidemic,ic,p)". How can I resolve it to get the plots?

My interest is to get the figure (3) i.e. The plot of nomalized sensitivity functions against time

function epidemic1()

clc
   
beta=0.4; gamma=0.04;
p1 = [beta gamma];
ic = [995;5;0];
F = ode45(@epidemic,ic,p);
    
   sol1 = solve(F,0,80)

   
figure(1)
        plot(sol1.Time,sol1.Solution,'-o')
        legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.4$, $\gamma=0.04$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')







Figure(2)
p2 = [0.2 0.1];
F.Parameters = p2;
F.Sensitivity = odeSensitivity;
sol2 = solve(F,0,80)

plot(sol2.Time,sol2.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')

U11 = squeeze(sol2.Sensitivity(1,1,:))'.*(p2(1)./sol2.Solution(1,:));
U12 = squeeze(sol2.Sensitivity(1,2,:))'.*(p2(2)./sol2.Solution(1,:));
U21 = squeeze(sol2.Sensitivity(2,1,:))'.*(p2(1)./sol2.Solution(2,:));
U22 = squeeze(sol2.Sensitivity(2,2,:))'.*(p2(2)./sol2.Solution(2,:));
U31 = squeeze(sol2.Sensitivity(3,1,:))'.*(p2(1)./sol2.Solution(3,:));
U32 = squeeze(sol2.Sensitivity(3,2,:))'.*(p2(2)./sol2.Solution(3,:));


Figure(3)
t = tiledlayout(3,2);
title(t,['Parameter Sensitivity by Equation','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel(t,'Time','Interpreter','latex')
ylabel(t,'\% Change in Eqn','Interpreter','latex')
nexttile
plot(sol2.Time,U11)
title('$\beta$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U12)
title('$\gamma$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U21)
title('$\beta$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U22)
title('$\gamma$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U31)
title('$\beta$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U32)
title('$\gamma$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
function dydt = epidemic(t,y,p)
   dydt = [0;0;0]; 
    S = y(1);
    I = y(2);
    R = y(3);
    N = S + I + R;

    dSdt = -beta*(I*S)/N;
    dIdt = beta*(I*S)/N - gamma*I;
    dRdt = gamma*I;

   
end
end

[Simon Diaz]

I found some issues. Here are the main :

• The ode solver should be called as oriented object fully, i think that help
• The ode function must have linked the parameter to gamma and beta
• The ode funcion must return the values of the slopes (it was fixed at dydt = [0;0;0])

Also figures must be called with lowercase (figure not Figure)

%% Corrected version epidemic1()

clc;close all;clear all
   
beta=0.4; gamma=0.04;
p1 = [beta gamma];
ic = [995;5;0];
    
fun_ode     = @(t,x,p)epidemic(t,x,p);
problem_ode = ode( ODEFcn = fun_ode );   % Define ODE function
problem_ode.InitialValue  = ic;            % Inivial value x(t=t0)
problem_ode.Parameters    = p1;

   sol1 = solve(problem_ode,0,80)

   
figure(1)
        plot(sol1.Time,sol1.Solution,'-o')
        legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.4$, $\gamma=0.04$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')


figure(2)
p2 = [0.2 0.1];
problem_ode.Parameters = p2;
problem_ode.Sensitivity = odeSensitivity;
sol2 = solve(problem_ode,0,80)

plot(sol2.Time,sol2.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')

U11 = squeeze(sol2.Sensitivity(1,1,:))'.*(p2(1)./sol2.Solution(1,:));
U12 = squeeze(sol2.Sensitivity(1,2,:))'.*(p2(2)./sol2.Solution(1,:));
U21 = squeeze(sol2.Sensitivity(2,1,:))'.*(p2(1)./sol2.Solution(2,:));
U22 = squeeze(sol2.Sensitivity(2,2,:))'.*(p2(2)./sol2.Solution(2,:));
U31 = squeeze(sol2.Sensitivity(3,1,:))'.*(p2(1)./sol2.Solution(3,:));
U32 = squeeze(sol2.Sensitivity(3,2,:))'.*(p2(2)./sol2.Solution(3,:));


figure(3)
t = tiledlayout(3,2);
title(t,['Parameter Sensitivity by Equation','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel(t,'Time','Interpreter','latex')
ylabel(t,'\% Change in Eqn','Interpreter','latex')
nexttile
plot(sol2.Time,U11)
title('$\beta$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U12)
title('$\gamma$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U21)
title('$\beta$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U22)
title('$\gamma$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U31)
title('$\beta$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U32)
title('$\gamma$, Eqn. 3','Interpreter','latex')
ylim([-5 5])


function dydt = epidemic(t,y,p)
 
    S = y(1);
    I = y(2);
    R = y(3);
    N = S + I + R;

    beta = p(1);
    gamma = p(2);

    dSdt = -beta*(I*S)/N;
    dIdt = beta*(I*S)/N - gamma*I;
    dRdt = gamma*I;

   dydt = [dSdt;dIdt;dRdt]; 
end

[SAHEED AJAO]

Thank you sir

Answer 1

[Star Strider]

I am not certain that this does everything you want, however it now has the virtue of running without error —

epidemic1

ans = 
  SolverID enumeration

    ode45

sol1 = 
  ODEResults with properties:

        Time: [0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80]
    Solution: [3x41 double]

sol2 = 
  ODEResults with properties:

           Time: [0 0.0027 27.4608 80]
       Solution: [3x4 double]
    Sensitivity: [3x2x4 double]

function epidemic1()

clc
   
beta=0.4; gamma=0.04;
p1 = [beta gamma];
ic = [995;5;0];
F = ode;
F.ODEFcn = @(t,y)epidemic(t,y,p1);
F.InitialValue = ic;
F.SelectedSolver
    
   sol1 = solve(F,0,80)

   
figure(1)
        plot(sol1.Time,sol1.Solution,'-o')
        legend('S','I','R')
title(["SIR Populations over Time","$\beta=0.4$, $\gamma=0.04$"],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')







figure(2)
p2 = [0.2 0.1];
F.Parameters = p2;
F.Sensitivity = odeSensitivity;
sol2 = solve(F,0,80)

plot(sol2.Time,sol2.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')

U11 = squeeze(sol2.Sensitivity(1,1,:))'.*(p2(1)./sol2.Solution(1,:));
U12 = squeeze(sol2.Sensitivity(1,2,:))'.*(p2(2)./sol2.Solution(1,:));
U21 = squeeze(sol2.Sensitivity(2,1,:))'.*(p2(1)./sol2.Solution(2,:));
U22 = squeeze(sol2.Sensitivity(2,2,:))'.*(p2(2)./sol2.Solution(2,:));
U31 = squeeze(sol2.Sensitivity(3,1,:))'.*(p2(1)./sol2.Solution(3,:));
U32 = squeeze(sol2.Sensitivity(3,2,:))'.*(p2(2)./sol2.Solution(3,:));


figure(3)
t = tiledlayout(3,2);
title(t,["Parameter Sensitivity by Equation","$\beta=0.2$, $\gamma=0.1$"],'Interpreter','latex')
xlabel(t,'Time','Interpreter','latex')
ylabel(t,'\% Change in Eqn','Interpreter','latex')
nexttile
plot(sol2.Time,U11)
title('$\beta$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U12)
title('$\gamma$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U21)
title('$\beta$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U22)
title('$\gamma$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U31)
title('$\beta$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U32)
title('$\gamma$, Eqn. 3','Interpreter','latex')
ylim([-5 5])

function dydt = epidemic(t,y,p)
   dydt = [0;0;0]; 
    S = y(1);
    I = y(2);
    R = y(3);
    N = S + I + R;

    dSdt = -beta*(I*S)/N;
    dIdt = beta*(I*S)/N - gamma*I;
    dRdt = gamma*I;

   
end
end

.

[SAHEED AJAO]

Thanks, I am grateful. I modified the code and ran it , I used Matlab R2013a to run it. I got this error: "Undefined function or variable 'ode'.

Error in epidemic11 (line 7)

F = ode; "

The modified code is given below

function epidemic11()
clc
   
beta=0.4; gamma=0.04;
p1 = [beta gamma];
ic = [995;5;0];
F = ode;
F.ODEFcn = @(t,y)epidemic(t,y,p1);
F.InitialValue = ic;
F.SelectedSolver
    
   sol1 = solve(F,0,80)
   
figure(1)
        plot(sol1.Time,sol1.Solution,'-o')
        legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.4$, $\gamma=0.04$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
figure(2)
p2 = [0.2 0.1];
F.Parameters = p2;
F.Sensitivity = odeSensitivity;
sol2 = solve(F,0,80)
plot(sol2.Time,sol2.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
U11 = squeeze(sol2.Sensitivity(1,1,:))'.*(p2(1)./sol2.Solution(1,:));
U12 = squeeze(sol2.Sensitivity(1,2,:))'.*(p2(2)./sol2.Solution(1,:));
U21 = squeeze(sol2.Sensitivity(2,1,:))'.*(p2(1)./sol2.Solution(2,:));
U22 = squeeze(sol2.Sensitivity(2,2,:))'.*(p2(2)./sol2.Solution(2,:));
U31 = squeeze(sol2.Sensitivity(3,1,:))'.*(p2(1)./sol2.Solution(3,:));
U32 = squeeze(sol2.Sensitivity(3,2,:))'.*(p2(2)./sol2.Solution(3,:));
figure(3)
t = tiledlayout(3,2);
title(t,['Parameter Sensitivity by Equation','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel(t,'Time','Interpreter','latex')
ylabel(t,'\% Change in Eqn','Interpreter','latex')
nexttile
plot(sol2.Time,U11)
title('$\beta$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U12)
title('$\gamma$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U21)
title('$\beta$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U22)
title('$\gamma$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U31)
title('$\beta$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U32)
title('$\gamma$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
function yprime = epidemic(~,y,~)
   yprime = [0;0;0]; 
    S = y(1);
    I = y(2);
    R = y(3);
    N = S + I + R;
    yprime(1) = -beta*(I*S)/N;
    yprime(2) = beta*(I*S)/N - gamma*I;
    yprime(3) = gamma*I;
end
end

[Sam Chak]

Aha, I see. The ode object was introduced in R2023b release, and the odeSensitivity feature is introduced in the latest release R2024a.

However, if you know the math equation (partial derivatives of the equations with respect to the parameters), then you can manually write a code for the Sensitivity function so that it runs in R2013a.

[Torsten]

The use of the "solve" command for an ODE structure is a very new MATLAB feature and not available for your old release.

[Star Strider]

@SAHEED AJAO — It would be best to upgrade. However, you can run it here (although not on MATLAB Online, since it reflects your version and installed Toolboxes).

As @Sam Chak mentioned, you can calculate it manually. If you have the Symbolic Math Toolbox, using the jacobian function and then the matlabFunction function makes that easier.

[Sam Chak]

The formula can be found here:

https://www.mathworks.com/help/matlab/math/sensitivity-analysis-of-epidemic-ode-parameters.html#EpidemicSAExample-5