%% Create a Rate Rule for the Rate of Change That Is Determined by Another Species
% This example shows how to create a rate rule where a species from one reaction can determine the rate of another reaction if it is in the second reaction rate equation. Similarly, a species from a reaction can determine the rate of another species if it is in the rate rule that defines that other species. Suppose you have a SimBiology model with three species (a, b, and c), one reaction (a -> b), and two parameters (k1 and k2). The rate equation is defined as b = -k_1 * a, and rate rule is dc/dt = k_2*a. The solution for the species in the reaction are:
% {\rm{a = }}{{\rm{a}}_{\rm{o}}}{{\rm{e}}^{{\rm{ - k_1t}}}}, {\rm{b = }}{{\rm{a}}_{\rm{o}}}{\rm{(1 - }}{{\rm{e}}^{{\rm{ - k_1t}}}}{\rm{)}}.
% Since the rate rule dc/dt = {k_2}*a is dependent on the reaction, dc/dt = {k_2}({a_o}{e^{ - {k_1}t}}). The solution is:
% c = {c_o} + {k_2}{a_o}/{k_1}(1 - {e^{ - {k_1}t}})
% Enter the following commands to set up a SimBiology model accordingly and simulate it.

m = sbiomodel('m');
c = addcompartment(m,'comp');
s1 = addspecies(m,'a','InitialAmount',10,'InitialAmountUnits','mole');
s2 = addspecies(m,'b','InitialAmount',0,'InitialAmountUnits','mole');
s3 = addspecies(m,'c','InitialAmount',5,'InitialAmountUnits','mole');
rxn = addreaction(m,'a -> b');
kl = addkineticlaw(rxn,'MassAction');
p1 = addparameter(kl,'k1','Value',1,'ValueUnits','1/second');
rule = addrule(m,'c = k2 * a','RuleType','rate');
kl.ParameterVariableNames = 'k1';
p2 = addparameter(m,'k2','Value',1,'ValueUnits','1/second');
[t,sd,species] = sbiosimulate(m);
plot(t,sd);
legend(species);
xlabel('Time');
ylabel('Species Amount');

% Copyright 2015 The MathWorks, Inc.