estimate the relative sensitivity coefficients

% TODO - Write the function declaration. Name the function SIRmodel
function dYdt = Lab3(t,Y,p) 
% TODO - Extract S1, S2
s1 = Y(1);
s2 = Y(2);
% TODO - Define the constants/parameters in the ODE system
V0 = p(1);
k1 = p(2);
V2 = p(3);
KM = p(4);
ds1dt = V0 - k1*s1;
ds2dt = k1*s1 - ((V2*s2)/(KM+s2));
% TODO - Create output column vector dYdt
dYdt = [ds1dt; ds2dt];
end

clear all
% Close all
% Timespan
tRange = [0 100];
% Initial condition
Y0 = [0;0];
% Parameters in ODE system
V0 = 1; % mM/min
k1 = 25; % 1/min
V2 = 2; % mM/min
KM = 5; % mM
% Want to choose a parameter p and a value between 1% and 5% so set vector
% to be P*[1;1.03]; so this give back the parameter the parameter with a 3% increase
delp = 0.01;
% For changing the param values in the loop
p_delp = V0*[1;1+delp];
% Run the model for the two values of the parameter p and p+delta_p
for i=1:length(p_delp)
    % Choose p and p + delta_p estimating the sensitivity
    V0 = p_delp(i);
    % Set the parameters to send to solver
    p = [V0,k1,V2,KM];
    % Call the ODE solver ode15s instead of ode45
    % to send parameters to the ode solver, use the following command:
    [tSol,YSol] = ode15s(@(tSol,YSol)Lab3(tSol,YSol,p),tRange,Y0);
    s1 = YSol(:,1);
    s2 = YSol(:,2);
    % Extract the steady state
     s1ss(i) = s1(end);
     s2ss(i) = s2(end);
end
% Compute dS1ss/dV0 and dS2ss/dV0
ds1ssdV0 = (s1ss(2) - s1ss(1))/delp;
ds2ssdV0 = (s2ss(2) - s2ss(1))/delp;
% Rescale the derivatives:
% new: (V0/s1ss)*ds1ss/dV0
% new: (V0/s2ss)*ds2ss/dV0
% Coefficients
scs1 = (V0./s1ss)*ds1ssdV0
scs2 = (V0./s2ss)*ds2ssdV0
% Plot solutions in time
figure(1)
clf
hold on
plot(tSol,YSol,'LineWidth',2)
xlabel('Time')
ylabel('Concentration')
legend('S_1','S_2')
set(gca,'FontSize',18)
grid on

Answer 1

[Torsten]

% TODO - Write the function declaration. Name the function SIRmodel
tRange = [0 100];
% Initial condition
Y0 = [0;0];
% Parameters in ODE system
V0 = 1; % mM/min
k1 = 25; % 1/min
V2 = 2; % mM/min
KM = 5; % mM
% Want to choose a parameter p and a value between 1% and 5% so set vector
% to be P*[1;1.03]; so this give back the parameter the parameter with a 3% increase
delp = 0.01;
% For changing the param values in the loop
p_delp = V0*[1;1+delp];
% Run the model for the two values of the parameter p and p+delta_p
for i=1:length(p_delp)
    % Choose p and p + delta_p estimating the sensitivity
    V0 = p_delp(i);
    % Set the parameters to send to solver
    p = [V0,k1,V2,KM];
    % Call the ODE solver ode15s instead of ode45
    % to send parameters to the ode solver, use the following command:
    [tSol,YSol] = ode15s(@(tSol,YSol)Lab3(tSol,YSol,p),tRange,Y0);
    % Extract the steady state
    s1ss(i) = YSol(end,1);
    s2ss(i) = YSol(end,2);
    % Save results
    t{i} = tSol;
    s1{i} = YSol(:,1);
    s2{i} = YSol(:,2);
end
% Compute dS1ss/dV0 and dS2ss/dV0
ds1ssdV0 = (s1ss(2) - s1ss(1))/delp;
ds2ssdV0 = (s2ss(2) - s2ss(1))/delp;
% Rescale the derivatives:
% new: (V0/s1ss)*ds1ss/dV0
% new: (V0/s2ss)*ds2ss/dV0
% Coefficients
scs1 = (p_delp(1)/s1ss(1))*ds1ssdV0
scs2 = (p_delp(1)/s2ss(1))*ds2ssdV0
% Plot solutions in time
figure(1)
plot(t{1},[s1{1},s2{1}],'LineWidth',2)
hold on
plot(t{2},[s1{2},s2{2}],'LineWidth',2)
xlabel('Time')
ylabel('Concentration')
legend('S_1','S_2')
set(gca,'FontSize',18)
grid on
end