Using simFunction Object in Simbiology when my stoichiometry depends on a model parameter

I have s SimBiology model in which the stochiometry of some equations depend on a model parameter called "DAR" as you can see below. I have DAR and some other paameters to change and see their effect on some observables. Once I change DAR and create a SimFunction, The stochiometry of equations do not change. Is there any way to use simFunction object and see the effect of DAR on the stchiometry. I can change the stochiometry in loop but I want it to be chnage as I use simFunction.

Thanks

reactionObj = addreaction(model,'Media.ADC_ext -> Media.Ab_ext + Media.PL_ext');

set (reactionObj, 'Stoichiometry', [-1 1 sbioselect(model,'Name','DAR').Value])

set (reactionObj, 'ReactionRate', 'kdec*Media.ADC_ext');

set (reactionObj, 'Notes', 'ADCs lose their payload, to produce unconjugated antibody and free payload');

set (reactionObj, 'Name', 'r1');

Answer 1

[Jeremy Huard]

Hi Mehdi,

the stoichiometry of reactions cannot be parametrized unfortunately.

But instead of using reactions you can use rate rules which will allow you to define stoichiometric coefficients as parameters.

Here is an example:

% model
model = sbiomodel("parstoich");
comp = addcompartment(model,"comp",1,Units="liter");

% species
addspecies(comp,"reactant",10,Units="milligram/liter");
addspecies(comp,"product",0,Units="milligram/liter");

% parameters including stoichiometric coeffs
n_r = addparameter(model,"n_r",3,Units="dimensionless");
addparameter(model,"n_p",2,Units="dimensionless");
addparameter(model,"k",1,Units="milligram/liter/hour");
addparameter(model,"mgperliter",1,Units="milligram/liter");

% ODEs
addrule(model,"reactant = -n_r*k*(reactant/mgperliter)^n_r","rate");
addrule(model,"product = n_p*k*(reactant/mgperliter)^n_r","rate");

% simulation with SimFunction
cs = getconfigset(model);
cs.TimeUnits = "hour";
cs.CompileOptions.UnitConversion = true;

equations = getequations(model)

equations = 
    'ODEs:
     d(reactant)/dt = -n_r*k*(reactant/mgperliter)^n_r
     d(product)/dt = n_p*k*(reactant/mgperliter)^n_r
     
     Parameter Values:
     n_r = 3 dimensionless
     n_p = 2 dimensionless
     k = 1 milligram/liter/hour
     mgperliter = 1 milligram/liter
     comp = 1 liter
     
     Initial Conditions:
     reactant = 10 milligram/liter
     product = 0 milligram/liter
     
     '

simfun = createSimFunction(model,"n_p",["reactant","product"],[],AutoAccelerate=false);

np_values = [2;4];
stopTime = 2;
sd = simfun(np_values,stopTime);

% plot
figure;
tl = tiledlayout('flow');
ax = gobjects(numel(np_values),1);
for i=1:numel(np_values)
    ax(i) = nexttile(tl);
    plot(ax(i),sd(i).Time,sd(i).Data);
    title(ax(i),"N_p = "+np_values(i))
end
lgd = legend(["reactant","product"]);
lgd.Layout.Tile = "east";
lgd.Box = "off";
xlabel(tl,"Time");
ylabel(tl,"mg/L");
set(ax,'XLimitMethod','padded','YLimitMethod','padded');
grid(ax,'on');
linkaxes(ax);

[Mehdi]

Hi Jeremy,

Thanks for the reply. This approach works for one reaction.

However, if we have a network of chemical reactions which have common species and in some reations, the stochiometry depends on some parameters, we should define rate rules for all those species. This could lead to mistakes. So, it would be better to use rate rules for all the parameters to avoid mistake. It seems now it is the only way.

Thanks

[Jeremy Huard]

I understand this can be tedious and error-prone for larger models.

But you can implement it fairly easily with the following procedure:

1. implement your model using reactions and some exemplary stoichiometric coefficients as you might have already done it
2. get the equations generated by SimBiology with getequations (or in the SimBiology Model Builder App)
3. create parameters for the stoichiometric coefficients
4. create rate rules by copying the ODE which need to be parametrized you got from #2
5. set the BoundaryCondition to true for the species that are now defined with rate rules

BoundaryCondition=true will have the species dynamics will be defined by the rate rule and not by the reactions.

You can set this property in Model Builder by selecting multiple species at once (by pressing Shift in Windows and clicking on the species) or in the command line with:

specs = sbioselect(model,"Type","species","Name",["reactant","product"]);

set(specs,"BoundaryCondition",true);

[Mehdi]

Thank you, make senses!

Answer 2

[Mehdi]

Hi,

I am working on a SimBiology code with Stochiometry depending on a parameter (DAR). The specise are calculated based on their concentrations (nmol/liter). The size of the component changes too which is based on repeated rules. I use SimFunction for paarmetrization, therefor as we discussed, I use rate rules for the specieses in which their concentration depends on DAR. The method was explained before in this thread. But I get the follwoing error. Is their any way to fix this error?

Rate rule for species 'PL_ext' is invalid. Species 'PL_ext' is in concentration and in a compartment with varying capacity. For this to be valid, the compartment capacity

must be varied by a rate rule.

[Arthur Goldsipe]

The error you report is intended to communicate that SimBiology does not allow the following the following combination of conditions:

• A species is in concentration.
• The concentration of that species is determined by a rate rule.
• The compartment (that the species is in) is determined by a repeated assignment rule or an algebraic rule.

One fix would be the one mentioned at the end of the error message: Replace the repeated assignment rule or algebraic rule with a rate rule that gives equivalent behavior.

Another option would be to change the species units to amount instead of concentration and use (species amount)/(compartment volume) anywhere you need to reference the concentration.

Also, one request: If you have additional questions, please post them as a new questions on MATLAB Answers. You posted this question in the "answer" field, and it makes it a bit harder to keep track of each question and the associated best answer.

[Arthur Goldsipe]

Oh, right after I provided the above answer I saw that you did repost this as a separate question. I'll repost my answer there.