## Definitions and Evaluations of Rules in SimBiology Models

### Overview

Rules are mathematical expressions that allow you to define or modify model quantities, namely compartment capacity, species amount, or parameter value.

Rules can take the form of initial assignments, assignments during the course of a simulation (repeated assignments), algebraic relationships, or differential equations (rate rules). Details of each type of rule are described next.

### Initial Assignment

An initial assignment rule lets you specify the initial value of a model quantity as a numeric value or as a function of other model quantities. It is evaluated once at the beginning of a simulation.

An initial assignment rule is expressed as ```
Variable =
Expression
```

, and the rule is specified as the
`Expression`

. For example, you could write an initial
assignment rule to set the initial amount of `species2`

to be
proportional to `species1`

as ```
species2 = k *
species1
```

where `k`

is a constant parameter.

### Repeated Assignment

A repeated assignment rule lets you specify the value of a quantity as a numeric value or as a function of other quantities repeatedly during the simulation. It is evaluated at every time step, which is determined by the solver during the simulation process.

A repeated assignment rule is expressed as ```
Variable =
Expression
```

, and the rule is specified as the
`Expression`

. For example, to repeatedly evaluate the total
species amount by summing up the species in different compartments, you could enter:
`xTotal = c1.X + c2.X`

, where `xTotal`

is a
nonconstant parameter, `c1`

and `c2`

are the name
of compartments where species x resides.

### Algebraic Rules

An algebraic rule lets you specify mathematical constraints on one or more model quantities that must hold during a simulation. It is evaluated continuously during a simulation.

An algebraic rule takes the form `0 = Expression`

, and the rule
is specified as the `Expression`

. For example, if you have a mass
conservation equation such as ```
species_total = species1 +
species2
```

, write the corresponding algebraic rule as ```
species1 +
species2 - species_total
```

.

However, repeated assignment rules are mathematically equivalent to algebraic rules, but result in exact solutions instead of approximated solutions. Therefore, it is recommended that you use repeated assignment rules instead of algebraic rules whenever possible. Use algebraic rules only when:

You cannot analytically solve the equations to get a closed-form solution. For example, there is no closed-form solution for

`x^4 + ax^3 + bx^2 + cx + k = 0`

whereas the closed-form solution for`kx – c = 0`

is`x = c/k`

.You have multiple equations with multiple unknowns, and they could be inconvenient to solve.

If you use an algebraic rule, rate rule, or repeated assignment to vary the value
of a parameter or compartment during the simulation, make sure the `ConstantValue`

property of the
parameter or `ConstantCapacity`

of the compartment
is set to `false`

.

### Repeated Assignment vs. Algebraic Rules

Repeated assignment rules are mathematically equivalent to algebraic rules, but result in exact solutions. However, algebraic rules are solved numerically, and the accuracy depends on the error tolerances specified in the simulation settings. In addition, there are several advantages to repeated assignment rules such as better computational performance, more accurate results since no rules have to be solved numerically (hence no approximations), and sensitivity analysis support.

**Tip**

If you can analytically solve for a variable, use a repeated assignment rule instead of an algebraic rule.

In repeated assignment rules, the constrained variable is explicitly defined as the left-hand side, whereas in algebraic rules it is inferred from the degrees of freedom in the system of equations. See also Considerations When Imposing Constraints.

### Rate Rules

A rate rule represents a differential equation and lets you specify the time derivative of a model quantity. It is evaluated continuously during a simulation.

A rate rule is represented as $$\frac{dVariable}{dt}=Expression$$, which is expressed in SimBiology as ```
Variable =
Expression
```

. For example, if you have a differential equation for
species *x*, $$\frac{dx}{dt}=k(y+z)$$, write the rate rule as: ```
x = k * (y +
z)
```

.

For more examples, see Rate Rule Examples.

### Evaluation Order of Rules

At the start of the simulation (that is, at simulation time = 0), SimBiology^{®} evaluates the initial assignment and repeated assignment rules as a
set of simultaneous constraints. SimBiology treats the rules as a unified system of
constraints and automatically reorders and evaluates them. The order in which the
rules appear in the model has no effect on the simulation results.

If a quantity is being modified by an assignment rule, the rule replaces initial
value properties, such as `InitialAmount`

,
`Capacity`

, or `Value`

. Similarly, a
variant altering such quantities has no effect because the value is superseded by
the assignment rules.

SimBiology throws an error if the model has circular dependencies in the initial assignment and repeated assignment rules. In other words, initial assignments and repeated assignments cannot have a variable that is explicitly or implicitly referenced on both the left- and right-hand sides of the equation.

For instance, you cannot create circular sets of assignments such as ```
a =
b + 1
```

and `b = a + 1`

, where *a* and
*b* are explicitly referenced on both sides of the equation. An
example of an implicit reference is when an assignment rule references a species in
concentration. In this case, the compartment that contains the species is implicitly
referenced.

**Warning**

You might observe different simulation results with respect to initial
assignments for previous releases of SimBiology (R2017a or earlier). To recover
the same simulation results at time = 0, as in R2017a or earlier, use the
`updateInitialAssignments`

function in the command line. If you
are using the **SimBiology** app, right-click the model from the
**Project Workspace** and select **Remove Order
Dependencies**.

### Conservation of Amounts During Simulation

During a simulation (that is, at simulation time > 0), SimBiology conserves species amounts rather than concentrations if there are any changes to the volume of a compartment where the species reside. In other words, if you have a repeated assignment rule or an event that changes the volume, then you see the effect of conservation of species amounts at time > 0.

However, at the beginning of a simulation (that is, at simulation time = 0), the concept of amount conservation does not apply because there are no changes before time = 0. Only one set of initial conditions exists and SimBiology uses the conditions at the start of the simulation. Specifically, at time = 0, SimBiology:

Initializes variables for species, compartments, and parameters using the corresponding

`InitialAmount`

,`Capacity`

, and`Value`

properties.Updates the values by replacing them with the corresponding alternate values from variants, if any.

Updates the values by evaluating initial assignment and repeated assignment rules as a set of simultaneous constraints. Therefore, the assignment rules replace initial values if model quantities are being modified by such rules or variants.

**Warning**

In previous releases (R2017a or earlier), if a repeated assignment changed a
compartment volume, SimBiology used the compartment capacity to determine the
initial amount and conserved it when the compartment volume changed at time = 0.
In R2017b or later, SimBiology uses the `InitialAmount`

property of the species as the initial condition at time = 0. Consider the
following model.

m = sbiomodel('m1') v = addcompartment(m,'v',10,'ConstantCapacity',0,'CapacityUnit','liter') p = addparameter(m,'p','ValueUnit','liter') r = addrule(m,'v = 100 * p','repeatedAssignment') s = addspecies(v,'s',50,'InitialAmountUnit','milligram/liter')

*s*as

```
50 milligram/liter * 10 liter = 500
milligram
```

, and then applied the repeated assignment rule
`v = 100 liter`

. So, the concentration of
*s*was then calculated and reported as

```
500
milligram/100 liter = 5 milligram/liter
```

at time = 0. In R2017b or later, SimBiology uses the `InitialAmount`

property of species *s*, and reports the initial amount of
*s* as `50 milligram/liter`

instead.

### Writing Rule Expressions

Use MATLAB^{®} syntax to write a mathematical expression for a rule. Note that no
semicolon or comma is needed at the end of a rule expression. If your algebraic,
repeated assignment, or rate rule expression is not continuous or differentiable,
see Using Events to Address Discontinuities in Rule and Reaction Rate Expressions before simulating your model.

### Considerations When Imposing Constraints

Suppose that you have a species `y`

whose amount is determined by
the equation `y = m * x - c`

. In SimBiology, the algebraic rule to
describe this constraint is written as `m * x - c - y`

. If you want
to use this rule to determine the value of `y`

, then
`m`

, `x`

, and `c`

must be
variables or constants whose values are known or determined by other equations.
Therefore, you must ensure that the system of equations is not overconstrained or
underconstrained. For instance, if you have more equations than unknowns, then the
system is overconstrained. Conversely, if you have more unknowns than the equations,
then the system is underconstrained.

**Tip**

The behavior of an underconstrained system could be fixed by adding additional
rules or by setting the `ConstantValue`

or `ConstantCapacity`

or `ConstantAmount`

property of some
of the components in the model.

### Rate Rule Examples

The following examples show how to create rate rules for different applications.

#### Create a Rate Rule for a Constant Rate of Change

This example shows how to increase the amount or concentration of a species by a constant value using the zero-order rate rule. For example, suppose species `x`

increases by a constant rate `k`

. The rate of change is:

$$dx/dt=k$$

Set the initial amount of species `x`

to 2, and the value of parameter `k`

to 1. Use the following commands to set up a SimBiology model accordingly and simulate it.

m = sbiomodel('m'); c = addcompartment(m,'comp'); s = addspecies(m,'x','InitialAmount',2); p = addparameter(m,'k','Value',1); r = addrule(m,'x = k','RuleType','rate'); [t,sd,species] = sbiosimulate(m); plot(t,sd); legend(species) xlabel('Time'); ylabel('Species Amount');

Alternatively, you could model a constant increase in a species using the Mass Action reaction `null`

-> `x`

with the forward rate constant `k`

.

clear m = sbiomodel('m'); c = addcompartment(m,'comp'); s = addspecies(m,'x','InitialAmount',2); r = addreaction(m,'null -> x'); kl = addkineticlaw(r,'MassAction'); p = addparameter(kl,'k','Value',1); kl.ParameterVariableNames = 'k'; [t,sd,species] = sbiosimulate(m); plot(t,sd); legend(species) xlabel('Time'); ylabel('Species Amount');

#### Create a Rate Rule for an Exponential Rate of Change

This example shows how to change the amount of a species similar to a first-order reaction using the first-order rate rule. For example, suppose the species `x`

decays exponentially. The rate of change of species `x`

is:

$dx/dt=-k*x$

The analytical solution is:

${C}_{t}={C}_{0}*{e}^{-kt}$

where ${C}_{t}$ is the amount of species at time t, and ${C}_{0}$ is the initial amount. Use the following commands to set up a SimBiology model accordingly and simulate it.

m = sbiomodel('m'); c = addcompartment(m,'comp'); s = addspecies(m,'x','InitialAmount',2); p = addparameter(m,'k','Value',1); r = addrule(m,'x = -k * x','RuleType','rate'); [t,sd,species] = sbiosimulate(m); plot(t,sd); legend(species); xlabel('Time'); ylabel('Species Amount');

If the amount of a species `x`

is determined by a rate rule and `x`

is also in a reaction, `x`

must have its `BoundaryCondition`

property set to `true`

. For example, with a reaction `a -> x`

and a rate rule $$\frac{dx}{dt}=k*x$$, set the `BoundaryCondition`

property of species `x`

to `true`

so that a differential rate term is not created from the reaction. The amount of `x`

is determined solely by a differential rate term from the rate rule. If the `BoundaryCondition`

property is set to `false`

, you will get the following error message such as `Invalid rule variable 'x' in rate rule or reaction`

.

#### Create a Rate Rule to Define a Differential Rate Equation

Many mathematical models in the literature are described with differential
rate equations for the species. You could manually convert the equations to
reactions, or you could enter the equations as rate rules. For example, you
could enter the following differential rate equation for a species
`C`

:

$$\frac{\text{dC}}{\text{dt}}\text{=vi-vdX}\frac{\text{C}}{\text{Kc+C}}\text{-kdC}$$

as a rate rule in SimBiology: `C = vi - (vd*X*C)/(Kc + C) - kd*C`

#### 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:

$$a={a}_{o}{e}^{-{k}_{1}t}$$, $$b={a}_{o}(1-{e}^{-{k}_{1}t})$$.

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');