Issue recreating result from a large system of ODEs

Question

Jorge Cossio 2019-4-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/457809-issue-recreating-result-from-a-large-system-of-odes

回答： darova 2019-4-23

I'm trying to recreate this figure from a modeling paper, and can't seem to get the same results.The figure I'm trying to recreate is the following (without the experimet data points).

The equations I pulled from the paper are the following. All constants were pulled from a table within the article, as well as the initial conditions. However, the results I get are strange.

Looking at the results from ode45 more than half of the values for the equations are NaN, and whatever numbers I do obtain are hundreds of orders of magnitude large which doesn't make sense. The plot I'm most interested is Cac, y(2).

Below is my code:

%V4 - rearranged equations, updated variables, all in units of M. Possibly
%need to separate into 2 systems and run y4-y7 with y1-y3 as inputs
clear all
close all
clc
%% Constants and other values
Rt = 4.4e4;
k1= 1.24; k2 = 2;k3 = 6.64; k4 = 5;k5 = 1; k6 = 100; k7 = 300;
kcicr = 1;
kcce = 0;
K1 = 0; K2 = 0.2;K3 = 0.15;K4 = 0.08;K5 = .321;
Khi = 0.38;
Kcicr = 0;
Bt = 120;
VcVs = 3.5;
Vp = 1.63;
Vex = 18.33;
Vhi = 4.76;
Qin_pas = 5;
Qin_stim = 2.5;
Ca0 = 0.1; %from paper parragraph
Ca_ex = 1000; %defined as 1mM whreas all other units are in uM, set to 1000uM
cs = 1; %***MISSING Cs INFO, have tried several values *********
%% Initial conditions for Cas and Cab
Ca_s0 = sqrt(k4/k5)*(Ca0/K3+Ca0);
Ca_b0 = (k6*Ca0/(k6*Ca0+k7))*Bt;
%% systems of eq
%for sys 5-6-7
kon = 5316.7e-6;
koff = 142.8;
kd = 0.12;
tau1 = 33;
tau2 = .005;
%% system 16-26-27-28 (MISSING qin)
f2 = @(t,y) [(-kon*y(1)*cs + koff* y(2));%eq 5 for Ru
            (kon*y(1)*cs - koff*y(2) - kd*y(2));%eq6 for Rb
            (kd*y(2));%eq 7 for R*
            (k1*y(3)*(y(5)/(K1+y(5)))-k2*y(4));%eq16
            (k3*((kcicr*y(5))/(Kcicr+y(5)))*(y(4)/(K2+y(4)))^3*y(6)-k4*(y(5)/K3+y(5))^2 +k5*y(6)^2 - k6*y(5)*(Bt - y(7)) + k7*y(7) - (Vp*y(5)/K4+y(5)) - (Vhi*y(5)/Khi*y(5)) - (Vex*y(5)/K5*y(5)) +(Qin_pas+Qin_stim));%eq 26 for Cac
            (VcVs*(kcce*(Ca_s0 - y(6))*(Ca_ex - y(6)) - (k3*(kcicr*y(5)/Kcicr+y(5))*((y(4)/K2+y(4))^3)*y(6) - k4*(y(5)/K3+y(5))^2 + k5*y(6)^2)));%eq27 for Cas
            (k6*y(5)*(Bt-y(7)) - k7*y(7))]; %eq28 for Cab
            
tspan = [0 700]; %timespan based on fig 5
init = [Rt 0 0 0 0 Ca_s0 Ca_b0]; %initial conditions for i, Cac, Cas, Cab,R
[t,y] = ode45(f2,tspan,init);
plot(t,y,'-o')
figure
%Cac plot
plot(t,y(:,5))

3 个评论
显示 1更早的评论隐藏 1更早的评论

Jorge Cossio 2019-4-22

编辑：Jorge Cossio 2019-4-22

Hi Star,

Thanks for your feedback. Indeed I went through and found several typos in the formulas I have gone back and corrected them. My results have sadly not changed.

I'm not sure what you mean regarding time derivatives such as

. Could you elaborate on the issue? In the model I'm following these were all used as constants, and their value was obtained from the paper's table of constants.

I've updated my code above, as well as labeled the equations in f2. Hopfully this addresses most of the issues you saw, and again thank you for your feedback.

RE: Edit

Star, I rearranged the system to follow the numerical order the equations are presented in their '()'. I rearranged the images in my original post to reflect this change as I could see how it was confusing. The system is still the same.

Star Strider 2019-4-22

My pleasure.

If the time derivatives are supposed to be constants, then use them as such. That was not obvious.

Note that NaN values are the result of 0/0, Inf/Inf, or any other NaN value already present. They will propagate through any recursive calculation (such as the numerical integration of differential equaitions), so the first NaN result causes all subsequent calculations using that value to be NaN as well.

See Debug a MATLAB Program (link) for help on detecting the problem.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

darova 2019-4-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/457809-issue-recreating-result-from-a-large-system-of-odes#answer_371824

在 MATLAB Online 中打开

Look at your (26) and (27). You missed parethensis and power

I'd suggest you to create separate function and write clearer code :

function du = myode(t,u)
% 
%   all constants
%
    % Exctract u() :
Ru = u(1);
Rb = u(2);
Rx = u(3);
i =  u(4);
Cac = u(5);
Cas = u(6);
Cab = u(7);
    % give names du() :
dRu = 
dRb = 
dRx = 
di = 
    % then equation (26) will be written:
dCac =  k3*( kcicr*Cac/(Kcicr+Cac) )*( i/(K2+i) )^3 *Cas ...
      - k4*( Cac/(K3+Cac) )^2 + k5*Cas^2 ...
      - k6 *Cac*(Bt-Cab) + k7*Cab ...
      - Vp *Cac^1.7/(K4^1.7 + Cac^1.7) ...             % parenthesis and power ^
      - Vhi*Cac^4.4/(Khi^4.4 + Cac^4.4) ...            % parenthesis and power ^
      - Vex*Cac /(K5 + Cac) + Qin_pas + Qin_stim;      % parenthesis
dCas = 
dCab = 
du = [dRu dRb dRx di dCac dCas dCab]';
end

Main code:

tspan = [0 700]; %timespan based on fig 5
init = [Rt 0 0 0 0 Ca_s0 Ca_b0]; %initial conditions for i, Cac, Cas, Cab,R
[t,y] = ode45(f2,tspan,init);
plot(t,y,'-o')
figure
%Cac plot
plot(t,y(:,5))

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Issue recreating result from a large system of ODEs

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Issue recreating result from a large system of ODEs

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论