MATLAB Answers

Find the steady states of a system.

26 views (last 30 days)
Ciaran
Ciaran on 7 Jan 2015
Answered: Ingrid Tigges on 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 Comments

Sign in to comment.

Accepted Answer

Ciaran
Ciaran on 7 Jan 2015
Basically the answer was to initialize the figure before the for loop

  0 Comments

Sign in to comment.

More Answers (1)

Ingrid Tigges
Ingrid Tigges on 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 Comments

Sign in to comment.

Communities

More Answers in the  SimBiology Community

Community Treasure Hunt

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

Start Hunting!

Translated by