Graph plotting of ODE with time variable

I would like to create a graph with time dependent, I have the below code:

S = 10e-6; %amount of streptavidin (μM)

B = 10e-6; %amount of biotin (μM)

S2 = S+B

S3 = S+2*B

S4 = S+3*B

S5 = S+4*B

syms S(t) S2(t) S3(t) S4(t) S5(t)

ode1 = diff(S,t) == (-4*kon*B*S)+koff*S2;

ode2 = diff(S2,t) == -(3*kon*B+koff)*S2+4*kon*B*S+2*koff*S3;

ode3 = diff(S3,t) == -(2*kon*B+2*koff)*S3+3*kon*B*S2+3*koff*S4;

ode4 = diff(S4,t) == -(kon*B+3*koff)*S4+2*kon*B*S3+4*koff*S5;

ode5 = diff(S5,t) == -4*koff*S5+kon*B*S4;

odes = [ode1; ode2; ode3; ode4; ode5]

Aodes = dsolve(odes)

where I tried to use the ode45 function but it seems doesn't work for me or maybe I misunderstand something. And the graph aims to show the value changes of S(t) to S5(t) from 0 to 200 mins.

The concept will be similar to below graph

[Chi Chun Jacky Yeung]

%% Change in free energy

% assume the temperature is 25 degree

koff = 2.7e-2; %1/s

kon = 1.9e4 % 1/Ms

Kd = kon/koff

Ka = 1/Kd

% Keq = Kd = 1/Ka

R = 8.314472; % J/K*mol universal gas constant

T = 25+273; % K

dG = log(Kd)*R*T % change in free energy

%% number of free and occupied binding site

S = 10e-6; %amount of streptavidin (μM)

B = 10e-6; %amount of biotin (μM)

S2 = S+B

S3 = S+2*B

S4 = S+3*B

S5 = S+4*B

syms S(t) S2(t) S3(t) S4(t) S5(t)

ode1 = diff(S,t) == (-4*kon*B*S)+koff*S2;

ode2 = diff(S2,t) == -(3*kon*B+koff)*S2+4*kon*B*S+2*koff*S3;

ode3 = diff(S3,t) == -(2*kon*B+2*koff)*S3+3*kon*B*S2+3*koff*S4;

ode4 = diff(S4,t) == -(kon*B+3*koff)*S4+2*kon*B*S3+4*koff*S5;

ode5 = diff(S5,t) == -4*koff*S5+kon*B*S4;

odes = [ode1; ode2; ode3; ode4; ode5]

Aodes = dsolve(odes)

Here's the full version of the code.

Answer 1

[Star Strider]

None of the initial conditions are defined (nor are ‘kon’ and ‘koff’), so the initial conditions will be replaced by symbolic constants that do not have numeric values. The result is that none of the equations can be evaluated with fplot or any other plot function.

[Chi Chun Jacky Yeung]

oh the kon and koff had been defined in the previous section, I just copied part of the code.

% assume the temperature is 25 degree

koff = 2.7e-2; %1/s

kon = 1.9e4 % 1/Ms

Kd = kon/koff

Ka = 1/Kd

% Keq = Kd = 1/Ka

R = 8.314472; % J/K*mol universal gas constant

T = 25+273; % K

dG = log(Kd)*R*T % change in free energy

%% number of free and occupied binding site

S = 10e-6; %amount of streptavidin (μM)

B = 10e-6; %amount of biotin (μM)

S2 = S+B

S3 = S+2*B

S4 = S+3*B

S5 = S+4*B

syms S(t) S2(t) S3(t) S4(t) S5(t)

ode1 = diff(S,t) == (-4*kon*B*S)+koff*S2;

ode2 = diff(S2,t) == -(3*kon*B+koff)*S2+4*kon*B*S+2*koff*S3;

ode3 = diff(S3,t) == -(2*kon*B+2*koff)*S3+3*kon*B*S2+3*koff*S4;

ode4 = diff(S4,t) == -(kon*B+3*koff)*S4+2*kon*B*S3+4*koff*S5;

ode5 = diff(S5,t) == -4*koff*S5+kon*B*S4;

odes = [ode1; ode2; ode3; ode4; ode5]

Aodes = dsolve(odes)

So here will be the full version

[Star Strider]

‘So here will be the full version’

Almost.

We still lack the initial conditions for them. They are a bit difficult to infer from the plot, since it is not obvious (at least to me) what those curves represent. There is only one in your system of differential equations, however all of the equations in the plot are labeled in terms of .

I find that confusing.

[Chi Chun Jacky Yeung]

I appologies that I didn't say it clear.

Basically the term refer to the S where S2 to S5 are refer to , , and resepctively. So after run the matlab it shows the result as below:

diff(S(t), t) == (27*S2(t))/1000 - (19*S(t))/25

diff(S2(t), t) == (19*S(t))/25 - (597*S2(t))/1000 + (27*S3(t))/500

diff(S3(t), t) == (57*S2(t))/100 - (217*S3(t))/500 + (81*S4(t))/1000

diff(S4(t), t) == (19*S3(t))/50 - (271*S4(t))/1000 + (27*S5(t))/250

diff(S5(t), t) == (19*S4(t))/100 - (27*S5(t))/250

which become '1x1symfunc', but when I use the ode45 with tspan[0 200] and y0 = 0 (as my understanding y0 is the first condition of the plot that I assume all of them start from 0), it doesn't work for me.

[Star Strider]

My best effort:

% assume the temperature is 25 degree
koff = 2.7e-2; %1/s
kon = 1.9e4; % 1/Ms
Kd = kon/koff;
Ka = 1/Kd;
% Keq = Kd = 1/Ka
R = 8.314472; % J/K*mol universal gas constant
T = 25+273; % K
dG = log(Kd)*R*T; % change in free energy
% %%  number of free and occupied binding site
S = 10e-6; %amount of streptavidin (μM)
B = 10e-6; %amount of biotin (μM)
S2 = S+B;
S3 = S+2*B;
S4 = S+3*B;
S5 = S+4*B;
S00 = 0.65;                                                         % Inferred From Plot
S20 = 0.08;                                                         % Inferred From Plot
S30 = 0.05;                                                         % Inferred From Plot
S40 = 0.05;                                                         % Inferred From Plot
S50 = 0.16;                                                         % Inferred From Plot
syms S(t) S2(t) S3(t) S4(t) S5(t)
ode1 = diff(S,t) == (-4*kon*B*S)+koff*S2;
ode2 = diff(S2,t) == -(3*kon*B+koff)*S2+4*kon*B*S+2*koff*S3;
ode3 = diff(S3,t) == -(2*kon*B+2*koff)*S3+3*kon*B*S2+3*koff*S4;
ode4 = diff(S4,t) == -(kon*B+3*koff)*S4+2*kon*B*S3+4*koff*S5;
ode5 = diff(S5,t) == -4*koff*S5+kon*B*S4;
odes = [ode1; ode2; ode3; ode4; ode5; S(0)==S00; S2(0)==S20; S3(0)==S30; S4(0)==S40; S5(0)==S50];
Aodes = dsolve(odes);
Ss = Aodes.S;
S2s = Aodes.S2;
S3s = Aodes.S3;
S4s = Aodes.S4;
S5s = Aodes.S5;

figure
fplot(Ss,[0 30])
hold on
fplot(S2s,[0 30])
fplot(S3s,[0 30])
fplot(S4s,[0 30])
fplot(S5s,[0 30])
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend('S','S_2','S_3','S_4','S_5', 'Location','NW')

There appear to be problems with the model, since the plot it produces does not match the plot you posted.

I must leave you to sort those, since I have no idea what you are doing.

[Chi Chun Jacky Yeung]

Thank you so much! it works!! And don't worry, because I use other sample of streptavidin and biotin to do the work so the result should be different.

[Star Strider]

As always, my pleasure!

[Chi Chun Jacky Yeung]

Btw I would like to know how you inferred the value? Cuz I changed some values so the value of S(0) to S5(0) will be different. Will you be able to show me how you did it? Thank you

[Star Strider]

I looked at the plot you posted, and estimated the initial conditions as the values at the beginning for each curve.

[Chi Chun Jacky Yeung]

Oooops, so is there no anyway to find out the value of S(0) to S5(0) by Matlab itself?

[Star Strider]

In certain situations, yes. If you have data for the variables and you want to fit them to a system of differential equations to estimate the parameters, you can certainly estimate the initial conditions as well. (See this Answer where I include them as the last elements of the parameter vector in order to estimate them along with the other parameters.)

Otherwise, MATLAB expects you to supply them, whether you are using the Symbolic Math Toolbox, or one of the numeric ODE solvers (such as ode45 or ode15s).