SimBiology

フォロー

I have a .sbml model stored in an .xml file. It contains a couple of observables. When I plto the simulation results, these are uniformly set to "1" when I use the ssa simulation, however, observables seems to work fine with ode simulations. This is the code: % Fetch model. model = sbmlimport('multistate.xml'); % Run simulation. cs = getconfigset(model); cs.SolverType = 'sundials'; cs.StopTime = 10.0; [t_ode, x_ode, names] = sbiosimulate(model); % Plot observables. plot(t_ode, x_ode(:,end-3:end)) xlabel('Time') ylabel('States') legend(names(end-3:end)) % Prepare model for SSA simulations (see https://uk.mathworks.com/matlabcentral/answers/1582584-simbio-running-stochastic-ssa-simulation-of-sbml-imported-model). parameters = sbioselect(model, 'Type', 'parameter'); for paramIndex = 1:numel(parameters) parameter = parameters(paramIndex); [~, usageTable] = findUsages(parameter); % Update any reaction that uses this parameter in its rate. if isempty(usageTable) % if check added by Torkel. continue end reactions = usageTable.Component(usageTable.Property == "ReactionRate"); for reactionIndex = 1:numel(reactions) reaction = reactions(reactionIndex); oldRate = reaction.ReactionRate; kineticLaw = reaction.KineticLaw; if isempty(kineticLaw) % Add a kinetic law kineticLaw = addkineticlaw(reaction, 'MassAction'); else % Update the existing kinetic law to mass action kineticLaw.KineticLawName = 'MassAction'; end kineticLaw.ParameterVariableNames = parameter.Name; newRate = reaction.ReactionRate; if ~strcmp(oldRate, newRate) warning("Reaction rate for reaction %s changed from %s to %s.", ... reaction.Name, oldRate, newRate); end end end % Run simulation. cs = getconfigset(model); cs.SolverType = 'ssa'; cs.StopTime = 10.0; [t_ssa, x_ssa, names] = sbiosimulate(model); % Plot observables. plot(t_ssa, x_ssa(:,end-3:end)) xlabel('Time') ylabel('States') legend(names(end-3:end)) % Non-observables plot still work. plot(t_ssa, x_ssa(:,1:end-4)) xlabel('Time') ylabel('States') legend(names(1:end-4)) % Observable plotting still work for the modified model when using ode. cs = getconfigset(model); cs.SolverType = 'sundials'; cs.StopTime = 10.0; [t_ode_2, x_ode_2, names] = sbiosimulate(model); plot(t_ode_2, x_ode_2(:,end-3:end)) xlabel('Time') ylabel('States') legend(names(end-3:end)) Which produces these plots: The observables can be plotted when the mode is simulated via ode: The observables are all "1" when simualted via ssa: The other species can be plotted without problem though: Finnaly, plotting the observables when the system has been simualted via ode still work, even after the model have been modified. There is something wrong where the observables values are not set when the ssa algorithm is used. Finally, I have attached the model file (extension changed to .txt so that mathworks website would accept it).

Hello, I have a biochemical reaction network model defined in the SBML format (attached). I can load the model using the sbmlimport: sbml_mod = sbmlimport('multistate.xml'); I can perform ODE simulations of this model using sbiosimulate: [time,x,names] = sbiosimulate(sbml_mod); However, if I want to perform stochastic SSA (GIllespie) simulations (https://uk.mathworks.com/help/simbio/ug/stochastic-solvers.html), I get problems. Running cs = getconfigset(sbml_mod,'active'); cs.SolverType = 'ssa'; cs.StopTime = 30; [t_ssa, x_ssa] = sbiosimulate(sbml_mod); I get errors: --> Error reported from Stochastic Compilation: Reaction named 'Reaction_2' has an empty kinetic law. For stochastic simulation all kinetic laws must be MassAction. (one for each of my 18 reactions) It turns out that when loading SBML models each reactions get no KineticLaw, which is needed for SSA simulations. E.g. sbml_mod.Reactions(1).KineticLaw simply returns []. However, all reactions in the SBML file are simple MassAction reactions, e.g. sbml_mod.Reactions(1).ReactionRate is 'kon*[R(a,l)]*[L(r)]'. I tried simply looping through all reactions and setting their KineticLaw field to 'MassAction'. But this still yields errors like --> Error reported from KineticLaw Validation: Parameter variable name on kinetic law '' is empty. The number of parameters on the kinetic law must match the number in its definition, and all parameter names must be set. So I somehow need to ftech the right parameters and att that to the MassAction kinetic law. It seems like this should be fairly straightforward (e.g. many other packages, like Copasi can do SSA simulations directly from SBML files). Is there some automatic way of setting all reactions to MassAction kinetic law, or some other way which would enable me to do an SSA simulation of the model?