Find the steady states of a system.

2 次查看（过去 30 天）
Ciaran
Ciaran 2015-1-7
回答： Ingrid Tigges 2015-1-7
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

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

采纳的回答

Ciaran
Ciaran 2015-1-7
Basically the answer was to initialize the figure before the for loop

更多回答（1 个）

Ingrid Tigges
Ingrid Tigges 2015-1-7
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.

社区

更多回答在  SimBiology Community

类别

在 Help Center 和 File Exchange 中查找有关 Common Operations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by