Matlab ode solver and Simbiology gives different results on the same model

Hello! I'm a PhD student from EPFL, I'm currently working with Matlab to try to model a simple double feedback loop in order to observe bistability, for this purpose I've built this set of ODEs:

function dy = model_equation(t,y)
dy = zeros(4,1);

sCDC2i=1.0;
dCDC2i=1.0;
sWEE1i=1.0;
dWEE1i=1.0;
a1=1.0;
a2=1.0;
b1=200.0;
b2=10.0;
v=1.0;
g1=4.0;
g2=4.0;
k1=30.0;
k2=1.0;

CDC2a=y(1);
CDC2i=y(2);
WEE1a=y(3);
WEE1i=y(4);

dy(1)=a1*CDC2i-b1*CDC2a*(v*WEE1a)^g1/(k1+(v*WEE1a)^g1);
dy(2)=sCDC2i-dCDC2i*CDC2i-a1*CDC2i+b1*CDC2a*(v*WEE1a)^g1/(k1+(v*WEE1a)^g1);
dy(3)=a2*WEE1i-b2*WEE1a*CDC2a^g2/(k2+CDC2a^g2);
dy(4)=sWEE1i-dWEE1i*WEE1i-a2*WEE1i+b2*WEE1a*CDC2a^g2/(k2+CDC2a^g2);

this set of ODEs is solved and plotted with the following commands:

function A = modelSim(time)
[T,Y]=ode15s(@model_equation,[0 time],[0 1.0 0 1.0]);
plot(T,Y(:,1),T,Y(:,2),T,Y(:,3),T,Y(:,4));hold on;
hleg1=legend('CDC2a','CDC2i','WEE1a','WEE1i');
set(hleg1,'Location','NorthEastOutside');
hold off;
A=[T,Y];

Now the problem is that the system is not bistable, but I've a species (CDC2a) that tends to accumulate and never reaches steady state. I was thinking that I had made some mistake while I was writing the model, but after several controls, I found that there aren't mistakes in the model. So I've tried to use simbiology, I copied and pasted the ODEs from the .m file to Simbiology, set all the parameters as in the ODEs file, and the model started to work correctly, I got exactly the results I wanted! So I started to check for difference between the Simbiology model and the .m file, this is what I've found:

• ODEs are exactly the same
• parameters and costants are exactly the same
• ode solver and its options are equal
• the ODEs system in the .m file is solved trough "[T,Y]=ode15s(@model_equation,[0 time],[0 1.0 0 1.0],options);" while in simbiology the given command is "data = sbiosimulate(m2, cs, [], []);" (where m2 is the model and cs the options for the solver)

Now I can't figure out what is happening, I've tried in many different ways to get those results to be the same but I wasn't able to do it. In order to continue my project I need to know at least why I obtain different results, and which one to trust. Any help will be much appreciated because at this point I'm running out of ideas! thank you very much in advance!

[Tolutola]

Have you tried the problem with simulink(using perhaps state space box or integrators)? It might shed more light on the problem.

[Jan]

Posting code as unordered list is very creative, but the standard formatting would be more convenient. Simply add an empty line before and after the code, mark the code and hit the "{} Code" button.

[Tiziano]

Although Im quite handy with Matlab and SimBiology, unfortunately I never learned to use simulink, so I don't know where to begin to simulate the model with this program, I tried to read the documentation but for one that has a biological background like mine, it seems a bit complicated. I will let you know if I will obtain some results from simulink.

[SARANYA SK]

how to solve simbiology ode equation in matlab

[SARANYA SK]

d/dt (EGFR) = 1/compartment*((compartment*(re1k1*[pro EGFR]-re1k2*[EGFR])) - (compartment*(re2 k1*[EGF]*[EGFR]-re2_k2*[L_EGFR])))

Answer 1

[Arthur Goldsipe]

Hi,

If you are using MATLAB version R2012b, you can get a text version of the ODE equations and see where they differ from your implementation. If your model is stored in a variable called model, just execute the command

equationText = getequations(model)

Since I don't have access to your SimBiology model or the code to create it, it's hard to say exactly how they might differ. However, if you look at the limiting behavior of the model you did provide, you can see why CDC2a keeps accumulating:

As time increases, CDC2a gets large, making the term CDC2a^g2/(k2+CDC2a^g2) approximately equal to 1. Likewise, as time increases, WEE1a is relatively small, making term (v*WEE1a)^g1/(k1+(v*WEE1a)^g1) approximately 0. Thus, your rate equations simplify to the following:

dy(1) = a*y2
dy(2) = sCDC2i - dCDC2i*CDC2i - a1*CDC2i
dy(3) = a2*WEE1i -b2*WEE1a
dy(4) = sWEE1i - dWEE1i*WEE1i -a2*WEE1i + b2*WEE1a

You can solve for a steady state for y(2:4) by solving for dy(2:4) = 0. Specificallly, at this steady state, y(2:4) = [.5, .1, 1]. However, at these conditions dy(1) = a*(0.1). Since this term is greater than zero, y(1) [or CMC2a] will continue to grow without bound, even though the other species remain constant.

-Arthur

[Tolutola]

Hi, Arthur is right. I created the same model in simulink and got the same results. Unfortunately, the "bad" results are correct.You need to be very careful when using tools like simbiology so as not to get misleading results.

[Arthur Goldsipe]

I also used the equations you posted to implement the model in SimBiology. I got the same results (CDC2a keeps increasing over time). So I think you must have made a mistake translating the SimBiology model to ODEs. I'm just brainstorming here, but did you leave off any stoichiometric coefficients? Did you correctly account for the effect of compartment volume? Did you convert all your units properly?

If it would be helpful, I can send you a copy of the SimBiology model I created. Just send me an email.

-Arthur