Perturbing Values and determing an SFO

I am using the following code

K1 = 0.005;

K2 = 0.005;

K3 = 0.005;

K4 = 0.005;

Kc = 0.5;

Kd = 0.02;

kdd = 0.01;

V2 = 1.5;

V4 = 0.5;

vd = 0.25;

VM1 = 3;

VM3 = 1.;

Synth = 0.025;

% Original Equations:

% PhsKnsKnt = [Synth-vd*X*(C/(Kd+C)) - kdd*C; M*VM3*((1-X)/(K3+(1-X))) - V4*

(X/(K4+X)); VM1*(C/(Kc+C))*((1-M)/(K1+(1-M))) - V2*(M/(K2+M))]; % Anonymous Function with substitutions:

PhsKnsKnt = @(T, P) [Synth - vd*P(2)*(P(1)/(Kd+P(1))) - kdd*P(1); P(3)*VM3*((1- P(2))/(K3+(1-P(2)))) - V4*(P(2)/(K4+P(2))); VM1*(P(1)/(Kc+P(1)))*((1-P(3))/(K1+(1- P(3)))) - V2*(P(3)/(K2+P(3)))];

[T P] = ode45(PhsKnsKnt, [0 100], [0.01; 0.01; 0.01]);

C = P(:,1);

X = P(:,2);

M = P(:,3);

figure(1)

plot(T, P)

legend('[C] (µ\itM\rm)', '[X] (µ\itM\rm)', '[M] (µ\itM\rm)', 'Location','NorthWest')

xlabel('Time (s)')

ylabel('Concentration')

grid

I am trying to perturb each of the values from K1- VM3 by 5%. For example for kd i tried to following : Kd = 0.02 -(0.02*0.05); I have changed the other constants accordingly by 5%. My issue is that my "C" value when running this code has yet to change even when i perturb different values. I am trying to figure out SFO to rank variables using the following:

%change in C / % change of the constant (being 5%).

Answer 1

[Star Strider]

Actually, it does change. It just doesn’t change very much.

Run this:

Tvct = linspace(0,33);
Kdv = [0.95  1  1.05]*Kd;
for k1 = 1:3
    Kd = Kdv(k1);
    PhsKnsKnt = @(T, P)  [Synth - vd*P(2)*(P(1)/(Kd+P(1))) - kdd*P(1);    P(3)*VM3*((1-P(2))/(K3+(1-P(2)))) - V4*(P(2)/(K4+P(2)));    VM1*(P(1)/(Kc+P(1)))*((1-P(3))/(K1+(1-P(3)))) - V2*(P(3)/(K2+P(3)))];
    [T P] = ode45(PhsKnsKnt, [Tvct],  [0.01;  0.01;  0.01]);
    CXM(:,:,k1) = P;
end

Csq = squeeze(CXM(:,1,:));
Csqr = [Csq(:,1)./Csq(:,2)    Csq(:,3)./Csq(:,2)];
Csqd = [Csq(:,1)-Csq(:,2)    Csq(:,3)-Csq(:,2)];

I put in Tvct to be sure the function was evaluated at the same time points (since the ode solvers are adaptive). I created Csq, a matrix of the three runs for C(t) with values for Kd varying from 0.95 to 1.05, and then calculated the ratios of them in Csqr. Csqd are the differences, since that may be what you want.

(I put up a transient earlier post about sensitivity functions. That works, but only on actual objective functions. I forgot for a minute that this is a DE. Doesn’t apply to DEs.)

[KayLynn]

not understanding the beginning of this...why does linspace go from 0,33 and kdv at 0.95 1 and 1.05....

Also in the ode solver line why is 0.01 listed three times....i meant to ask that originally

[Star Strider]

I restricted the time to one oscillation here because it’s easier to see the details that way. Kdv of 0.95*Kd = Kd-(Kd*0.5), similarly for 1.05.

The initial conditions of the DE in your original post were all 0.01. They have to be specified for each variable in the MATLAB ode solvers, thus the [3 x 1] vector.

P.S. — Change the solver from ode45 to ode15s. It’s a bit better at these sorts of DEs.

P.P.S. — Plot this to see the effect:

figure(2)
plot(T,Csqr)
legend('[C] (µ\itM\rm) K_d = 0.95xK_d',  '[C] (µ\itM\rm) K_d = 1.05xK_d', 'Location','NorthWest')
grid

Answer 2

[KayLynn]

I do not have access to the tool box. My issue seems to be that my "c" varaible only changes when I add the perturbation rather than subtract..im not sure if that helps or changes things.