Find the steady states of a system.
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Ciaran
am 7 Jan. 2015
Beantwortet: Ingrid Tigges
am 7 Jan. 2015
I am building a SimBiology Model. I want to simulate from randomly generated initial conditions numerous times and find the steady state of a particular specie (PLS). MY current code does this in a graphical sense but I need the actual numbers in a vector (for example) for subsequent use. Does anybody know how to do this? My current code is:
%Random Number Generations and Assignment
No_of_simulations = 10; %i.e. points
number_of_parameters = 4; %i.e. dimenions
lb=[1e-3];
ub=[0.1];
x = lhsdesign(number_of_parameters ,No_of_simulations);
D = bsxfun(@plus,lb,bsxfun(@times,x,(ub-lb)));
for i=x,
n1 = i(1);
n2 = i(2);
n3 = i(3);
n4 = i(4);
%Model generation.
Mobj = sbiomodel('u-PA_parameter_generation');
comp_obj = addcompartment(Mobj, 'plasma');
%Reaction 1
Robj1 = addreaction(Mobj, 'Pro-u-PA + PLG -> PLS + Pro-u-PA');
Kobj1= addkineticlaw(Robj1, 'MassAction');
Pobj1 = addparameter(Kobj1, 'keff_zymogen',0.035);
set(Kobj1, 'ParameterVariableNames','keff_zymogen');
%Reaction 2
Robj2 = addreaction(Mobj, 'PLS + Pro-u-PA -> PLS + u-PA');
Kobj2= addkineticlaw(Robj2, 'MassAction');
Pobj2 = addparameter(Kobj2, 'keff_PLS',40);
set(Kobj2, 'ParameterVariableNames','keff_PLS');
%Reaction 3
Robj3 = addreaction(Mobj, 'u-PA + PLG -> u-PA + PLS');
Kobj3= addkineticlaw(Robj3, 'MassAction');
Pobj3 = addparameter(Kobj3, 'keff_pos',0.9);
set(Kobj3, 'ParameterVariableNames','keff_pos');
%Reaction 4
Robj4 = addreaction(Mobj, 'Pro-u-PA -> null');
Kobj4= addkineticlaw(Robj4, 'MassAction');
Pobj4 = addparameter(Kobj4, 'u1',0.084);
set(Kobj4, 'ParameterVariableNames','u1');
%Reaction 5
Robj5 = addreaction(Mobj, 'PLG -> null');
Kobj5= addkineticlaw(Robj5, 'MassAction');
Pobj5 = addparameter(Kobj5, 'u2',0.032);
set(Kobj5, 'ParameterVariableNames','u2');
%Reaction 6
Robj6 = addreaction(Mobj, 'PLS -> null');
Kobj6= addkineticlaw(Robj6, 'MassAction');
Pobj6 = addparameter(Kobj6, 'u1',0.084); %Same as R4
set(Kobj6, 'ParameterVariableNames','u1');
%Reaction 7
Robj7 = addreaction(Mobj, 'u-PA -> null');
Kobj7= addkineticlaw(Robj7, 'MassAction');
Pobj7 = addparameter(Kobj7, 'u1',0.084); %Same as R4
set(Kobj7, 'ParameterVariableNames','u1');
%Reaction 8
Robj8 = addreaction(Mobj, 'null -> Pro-u-PA');
Kobj8= addkineticlaw(Robj8, 'MassAction');
Pobj8 = addparameter(Kobj8, 'a1',0.0032);
set(Kobj8, 'ParameterVariableNames','a1');
%Reaction 9
Robj9 = addreaction(Mobj, 'null -> PLG');
Kobj9= addkineticlaw(Robj9, 'MassAction');
Pobj9 = addparameter(Kobj9, 'a2',0.01);
set(Kobj9, 'ParameterVariableNames','a2');
%setting species concentrations
Sobj1 = sbioselect(Mobj,'Type','species','Name','Pro-u-PA');
set(Sobj1, 'InitialAmount',n1)
Sobj2 = sbioselect(Mobj,'Type','species','Name','PLG');
set(Sobj2, 'InitialAmount',n2)
Sobj3 = sbioselect(Mobj,'Type','species','Name','PLS');
set(Sobj3, 'InitialAmount',n3)
Sobj4 = sbioselect(Mobj,'Type','species','Name','u-PA');
set(Sobj4, 'InitialAmount',n4)
%simulate
config = getconfigset(Mobj);
set(config,'StopTime',1000)
[t_ode, x_ode, names] = sbiosimulate(Mobj);
hold on
figure;
set(gcf, 'color','white');
plot(t_ode, x_ode(:,1:end));
legend(names)
end
Thanks in Advance
0 Kommentare
Akzeptierte Antwort
Weitere Antworten (1)
Ingrid Tigges
am 7 Jan. 2015
Given that the steady state is defined as dx/dt=0 which is approximately delta x/delta t you can check whether the difference in x of two consecutive time points is 0
diff_x= diff(x_ode);
Due to numeric errors you want to have those differences that are below a certain threshold:
threshold = 1e-15;
x_equals_0 = diff_x<threshold;
The indices of those values in x_equals_0 that are 1 belong to points where dx/dt=0.
0 Kommentare
Communitys
Weitere Antworten in SimBiology Community
Siehe auch
Kategorien
Mehr zu Extend Modeling Environment finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!