can I directly change the ODE?

I want to ensure that I truly understand how to use Simbiology by reconstructing an old model. I have the diagram, parameters (rates), initial values, ODEs, and result plots. After using this information to rebuild the model, I found that my results differ from the original ones. I would like to know if any other factors can affect the simulation results besides the data mentioned above. Can I directly modify the ODEs?
Thank you!

8 Kommentare

Walter Roberson
Walter Roberson am 20 Jul. 2023
Is there a choice of ode solver?
Lo
Lo am 20 Jul. 2023
I use ode15s
Sam Chak
Sam Chak am 20 Jul. 2023
Hi @Lo,
Can you show (i) the ODEs as well as (ii) the expected results/plots from original model based on the same initial values?
Lo
Lo am 20 Jul. 2023
they ode
my ode
they result
my result (StopTime : 60 mins, unitconversion : false)
my result (StopTime : 60 mins, unitconversion : ture, DefaultSpeciesDimension : concentration)
Sam Chak
Sam Chak am 20 Jul. 2023
Hi @Lo, Can you paste the copyable code here?
Also provide the initial values of the 12 states and constants of the parameters from k_padlock to k8.
Lo
Lo am 20 Jul. 2023
sorry I don't know how view my model in code form, but I can first provide the parameter and value
parameter
initial value
Thanks @Lo. Only the ODEs are left.
Can you type out or paste the 12 ODEs like this? It is for testing the ODEs in MATLAB.
dM/dt = - k_payload*M*P - Mdeg*M
dP/dt = - k_payload*M*P
...
Lo
Lo am 21 Jul. 2023
like this ways ?
d(M)/dt = -(k_padlock*M*P) - (Mdeg*M)
d(P)/dt = -(k_padlock*M*P)
d(MP)/dt = (k_padlock*M*P) - (Vmcircular*MP/(Kmcircular+MP))
d(SplR)/dt = -(Vmcircular*MP/(Kmcircular+MP))
d([MP-circular])/dt = (Vmcircular*MP/(Kmcircular+MP)) - (Vmphi29*[MP-circular]/(Kmphi29+[MP-circular]))
d(phi29)/dt = -(Vmphi29*[MP-circular]/(Kmphi29+[MP-circular]))
d([R-free])/dt = (Vmphi29*[MP-circular]/(Kmphi29+[MP-circular])) - (RCPdeg*[R-free]) - (k1*[R-free]*[F-free]-k2*RF) - (k3*[R-free]*[Q-free]-k4*RQ)
d([F-free])/dt = -(k1*[R-free]*[F-free]-k2*RF) - (k7*[F-free]*RQ-k8*RFQ)
d(RF)/dt = (k1*[R-free]*[F-free]-k2*RF) - (k5*RF*[Q-free]-k6*RFQ)
d([Q-free])/dt = -(k3*[R-free]*[Q-free]-k4*RQ) - (k5*RF*[Q-free]-k6*RFQ)
d(RQ)/dt = (k3*[R-free]*[Q-free]-k4*RQ) - (k7*[F-free]*RQ-k8*RFQ)
d(RFQ)/dt = (k5*RF*[Q-free]-k6*RFQ) + (k7*[F-free]*RQ-k8*RFQ)

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Fulden Buyukozturk
Fulden Buyukozturk am 20 Jul. 2023
Bearbeitet: Walter Roberson am 20 Jul. 2023

1 Stimme

You cannot edit ODEs directly but you can inspect them using getequations from the command-line or Show model equations option from Model Builder. You can also inspect quantity initial conditions by selecting Initial Conditions option on the Home tab.
Also, please see this page that describes how SimBiology derives ODEs: https://www.mathworks.com/help/simbio/ug/deriving-simbiology-odes-equations-from-reactions.html
Sam Chak
Sam Chak am 21 Jul. 2023

0 Stimmen

Hi @Lo
The initial values for MP0 and MPcircular0 are not provided. So I made up some values and it displays the trend similar to the following:
tspan = [0 600];
M0 = 0.5;
P0 = 0.05;
MP0 = 0.5; % not provided
SplR0 = 25;
MPcircular0 = 0; % not provided
phi290 = 12.5;
Rfree0 = 0;
Ffree0 = 180664245;
RF0 = 0;
Qfree0 = 361328490;
RQ0 = 0;
RFQ0 = 0;
x0 = [M0 P0 MP0 SplR0 MPcircular0 phi290 Rfree0 Ffree0 RF0 Qfree0 RQ0 RFQ0];
% options = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);
[t, x] = ode15s(@odefcn, tspan, x0);
% Solution Plot
plot(t/60, x(:,12), 'linewidth', 1.5),
grid on, xlabel('Time / minutes'), ylabel('RFQ')
function xdot = odefcn(t, x)
xdot = zeros(12, 1);
% definitions
M = x(1);
P = x(2);
MP = x(3);
SplR = x(4);
MPcircular = x(5);
phi29 = x(6);
Rfree = x(7);
Ffree = x(8);
RF = x(9);
Qfree = x(10);
RQ = x(11);
RFQ = x(12);
% parameters
k_padlock = 32000;
Mdeg = 0.000019254;
Vmcircular = 0.008;
Kmcircular = 1;
Vmphi29 = 112;
Kmphi29 = 31;
RCPdeg = 0;
k1 = 1e12;
k2 = 0;
k3 = 1e12;
k4 = 0;
k5 = 1e12;
k6 = 0;
k7 = 1e12;
k8 = 0;
% dynamics
xdot(1) = - k_padlock*M*P - Mdeg*M;
xdot(2) = - k_padlock*M*P;
xdot(3) = k_padlock*M*P - (Vmcircular*MP)/(Kmcircular + MP);
xdot(4) = - (Vmcircular*MP)/(Kmcircular + MP);
xdot(5) = (Vmcircular*MP)/(Kmcircular + MP) - (Vmphi29*MPcircular)/(Kmphi29 + MPcircular);
xdot(6) = - (Vmphi29*MPcircular)/(Kmphi29 + MPcircular);
xdot(7) = (Vmphi29*MPcircular)/(Kmphi29 + MPcircular) - RCPdeg*Rfree - (k1*Rfree*Ffree - k2*RF) - (k3*Rfree*Qfree - k4*RQ);
xdot(8) = - (k1*Rfree*Ffree - k2*RF) - (k7*Ffree*RQ - k8*RFQ);
xdot(9) = (k1*Rfree*Ffree - k2*RF) - (k5*RF*Qfree - k6*RFQ);
xdot(10) = - (k3*Rfree*Qfree - k4*RQ) - (k5*RF*Qfree - k6*RFQ);
xdot(11) = (k3*Rfree*Qfree - k4*RQ) - (k7*Ffree*RQ - k8*RFQ);
xdot(12) = (k5*RF*Qfree - k6*RFQ) + (k7*Ffree*RQ - k8*RFQ);
end

5 Kommentare

Lo
Lo am 21 Jul. 2023
I notice that the number of RFQ seems not right or I just misunderstand the figure ?
Sam Chak
Sam Chak am 21 Jul. 2023
Bearbeitet: Sam Chak am 21 Jul. 2023
Hi @Lo
I do not understand about the physical meaning of this 12th-order Simbiology system. The result is directly obtained from mathematically solving the ODEs based on the initial values and the parameters that you supplied, so that you can compare it with the result from the Simbiology Toolbox.
By the way, what is the right number of RFQ?
Lo
Lo am 23 Jul. 2023
the right number of RFQ should be near 10^-4, maybe 9.xxx*10^-4.
Sam Chak
Sam Chak am 23 Jul. 2023
Hi @Lo
The final value of RFQ in the simulation depends on the initial values of all 12 states. But you only gave 10 initial values (scoll up). The other two initial values of MP0 and MPcircular were not provided. What are they?
Please insert the values in the MATLAB code and run it. Then show your result of the RFQ here.
Lo
Lo am 2 Aug. 2023
I apologize for the late response due to being busy. The result I obtained is exactly the same as yours, but I still don't know what caused the difference in results before.

Melden Sie sich an, um zu kommentieren.

Communitys

Weitere Antworten in  SimBiology Community

Kategorien

Tags

Gefragt:

Lo
am 20 Jul. 2023

Kommentiert:

Lo
am 2 Aug. 2023

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by