How should I modify the rate rule using Symbiology so that it works?

1 view (last 30 days)
I am attempting to model a system for the Guerbert reaction. The reaction begins with Ethanol, which couples with itself to form Butanol, which couples with Ethanol to form Hexanol....and so on.
Here's the code:
m1 = sbiomodel('Guerbert Reaction');
c1 = addcompartment(m1, 'Plug-flow reactor');
% compartment added
s1 = addspecies(m1, 'EtOH');
s2 = addspecies(m1, 'BtOH');
s3 = addspecies(m1, 'HxOH');
% create all species
set(s1, 'InitialAmount', 5)
set(s2, 'InitialAmount', 0)
set(s3, 'InitialAmount', 0)
%set initial amounts for all species
r1 = addreaction(m1, 'EtOH + EtOH -> BtOH');
r2 = addreaction(m1, 'EtOH + BtOH -> HxOH');
%reactions added, parent object is m1 model
p1 = addparameter(m1, 'k2_2', 'Value', 0.1);
p2 = addparameter(m1, 'k2_4', 'Value', 0.4);
%parameters added to model object. This ensures that all objects within the model can use it
k1 = addkineticlaw(r1, 'MassAction');
k2 = addkineticlaw(r2, 'MassAction');
% kinetic laws added, of Mass Action type. Parent is the reaction object
p1 = addparameter(k1, 'k2_2', 'Value', 0.1);
p2 = addparameter(k2, 'k2_4', 'Value', 0.4);
%parameters added for use within the kinetic law
setparameter(k1, 'Forward Rate Parameter', 'k2_2')
setparameter(k2, 'Forward Rate Parameter', 'k2_4')
% defining the kinetic law by specifying which parameters will be substituted in the expression
ru1 = addrule(m1, 'EtOH = (-k2_2 * EtOH)/(1+EtOH) - (k2_4 * BtOH)/(1+EtOH)')
ru2 = addrule(m1, 'BtOH = (k2_2 * EtOH)/(1+EtOH) - (k2_4*BtOH)/(1+EtOH)')
ru3 = addrule(m1, 'HxOH = (k2_4 * BtOH)/(1*EtOH)')
% rate rules
set(ru1, 'RuleType', 'Rate')
set(ru2, 'RuleType', 'Rate')
set(ru3, 'RuleType', 'Rate')
% rule type is "rate" because the rules are differential}
However, when I verify the model, I get this error message:
>> verify(m1)
Error using SimBiology.Model/verify
--> Error reported from ODE Compilation:
Invalid rule variable 'EtOH' in rate rule 'EtOH = (-k2_2 * EtOH)/(1+EtOH) - (k2_4 * BtOH)/(1+EtOH)'. This object is being used in a reaction or in another rule.
What does this mean? How may I correct this?

Answers (1)

Jeremy Huard
Jeremy Huard on 30 May 2022
This is an old question but an answer could still benefit others who bumped into the same issue.
The key piece of the error message is: This object is being used in a reaction or in another rule
Indeed, the ODE driving the dynamics of a species can only be defined by EITHER a set of reactions OR a rate rule. But not both.
For the code above to work, 3 options can be considered:
  1. remove the reactions from the model and only keep the rate rules
  2. remove the rate rule and only keep the reactions. There one can use a built-in kinetic law or type in a custom reaction rate.
  3. one could keep the reactions for visualization purposes and the rate rules. but for this to work, the species property BoundaryCondition must be set to true.
Side note:
There is no need to add the parameters to the kinetic laws (lines of code with addparameter(kx,...) ). The forward rate parameters can be the model-scoped parameters.
The code above creates two independent k2_2 parameters: one scoped to the model (accessible by all reactions), one scoped to the reaction (only accessible by this reaction). Changing one does not change the other.
Here, both lines of code with addparameter(kx,...) can be removed.

Community Treasure Hunt

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

Start Hunting!

Translated by