Issue recreating result from a large system of ODEs

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))

[Star Strider]

It appears to me that you are using and (and the others) interchangably in your code, even though they have different values. You need to go through it and proofread it carefully. I found several instances in which your code differed from the images of your equations in that and other respects. Several expressions that should be raised to powers are not, and time derivatives and others are constants.

The errors are too numerous for me to trace them and correct them.

EDIT —

You edited your code to be an entirely different system of equations. They do not appear to be the same as the equations in the image.

[Jorge Cossio]

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]

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]

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))