MATLAB: Coupled rate ODEs with ode45

differential equationskinetic modelsode23ode45rate equations

Dear all,

I have been having trouble trying to model the concentrations of 8 species in the attached file showing the reactions. I have the rate constant values, initial conditions and a time frame but for some reason the plot seems to display no changes in concentration. I have attached my function and script. Any pointers or tips would be greatly appreciated.

Thanks in advance.

Best Answer

The concentrations change appropriately, however they don’t change much and the concentrations are vanishingly small. That’s the reason they don’ appear to change.

Increasing the final time (by 15 times) demonstrates saturation kinetics:

[t,state]=ode15s(@reactions,[0 51.5*15],[2.71E-5;7.8E-6;9.2E-6;0;1.0614E-4;3.359E-3;3.9278E-4;4.521E-4]);
reactants = {'TAG','DAG','MAG','GLY','FAEE','ET','FFA','W'};
figure
for k = 1:8
    subplot(4,2,k)
    plot(t,state(:,k))
    xlabel('Time')
    ylabel('Concentration')
    title(sprintf('$C_{%s}(t)$',reactants{k}), 'Interpreter','latex')
    grid
end

producing:

Adding:

pos = get(gcf,'Position')
set(gcf, 'Position',pos+[0 -500 200 500])

after the loop enlarges the figure making the subplot plots easier to see. (I did not use that in the posted figure.)

Related Solutions

MATLAB: Modelling a system of ODEs and need to add inputs after t=0

I am not exactly certain what you are doing or what your pharmacokinetic model is. Not knowing more, I would do something like this (assuming ‘t’ is the time vector and the integrated values are stored in matrix ‘y’):

Use first set of initial conditions ‘ic’

tspan = linspace(0, 10, 10);

Run ode15s with your ODE function

ic = y(end,:);
t0 = t(end);
ic(appropriate_compartment) = dose(1);

Use new set of initial conditions ‘ic’

tspan = linspace(t0, t0+60, 60);

Run ode15s with your ODE function

ic = y(end,:);
t0 = t(end);
ic(appropriate_compartment) = dose(2);

... continuing for as many steps as necessary. Adjust the number of elements in ‘tspan’ with the third argument of linspace. You will need to vertically concatenate the ‘t’ vectors and ‘y’ matrices from each step (loop iteration) to create the complete simulation.

I’m simply sketching this, corresponding to similar problems I’ve worked on.

Apparently the SimBiology has a set of functions that will make alll this straightforward, although I have no experience with it. (I don’t do such simulations frequently enought to justify buying it.)

MATLAB: Finding min and max of ODEs

See this example

% defining the time range:
T=12;
trange=[0 T];
% defining initial values (concentrations) for c,s,e i p:
cs0 = [0; 1; 1; 0];
  % [t,x] = ode45(@reaction, trange, initial values)
[t,unknown]=ode45(@reaction,trange,cs0);
c = unknown(:, 1);
s = unknown(:, 2);
e = unknown(:, 3);
p = unknown(:, 4);
[c_max, idx_c_max] = max(c);
t_c_max = t(idx_c_max);
[e_min, idx_e_min] = min(e);
t_e_min = t(idx_e_min);
% plot drawing
hold on
plot(t,c,'m',t,s,'c',t,e,'r',t,p,'k')
plot(t_c_max, c_max, 'm+', t_e_min, e_min, 'r+');
ylim([0 1.1])
xlabel('time, t')
ylabel('concentration')
legend({' enzyme-substrate complex concentration',' substrate concentration',' enzyme concentration ',' product concentration'}, "Location","best")
function dydt = reaction(~, cse)
% Defining reaction parameters:
k1 = 6;
k2 = 1;
kminus1 = 2;
% Defining reactions
reaction=zeros(4,1);
reaction(1) = k1*cse(2)*cse(3)-kminus1*cse(1)-k2*cse(1);
reaction(2) = -k1*cse(2)*cse(3)+kminus1*cse(1);
reaction(3) = -k1*cse(2)*cse(3)+kminus1*cse(1)+k2*cse(1);
reaction(4) = k2*cse(1);
dydt=[reaction(1); reaction(2); reaction(3); reaction(4)];
end

Best Answer

Related Solutions

MATLAB: Modelling a system of ODEs and need to add inputs after t=0

MATLAB: Finding min and max of ODEs

Related Question