Find the steady states of a system.
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Plus de réponses (1)
Ingrid Tigges
le 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 commentaires
Communautés
Plus de réponses dans SimBiology Community
Voir également
Catégories
En savoir plus sur Extend Modeling Environment dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!