Find the steady states of a system.

5 次查看(过去 30 天)
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 CenterFile Exchange 中查找有关 Extend Modeling Environment 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by