Use of Basic Elements II – Ordinary Differential Equation System

On this page

This second example shall illustrate the formulation of an ordinary differential equation system (ODE system) in MOSAICmodeling. This example is the “Oregonator”, a theoretical model for a type of autocatalytic reactions. In this example, we will use the publication by D’Ambrosio et al. (2018), DOI: 10.1007/s10910-018-0922-5, as basis.

Model description

The model consists of three differential equations:

$\begin{align*}&\frac{\mathrm{d}X}{\mathrm{d}t} = k_1 \cdot A \cdot Y - k_2 \cdot X \cdot Y + k_3 \cdot A \cdot X - 2 \cdot k_4 \cdot X^2 \\[2ex]&\frac{\mathrm{d}Y}{\mathrm{d}t} = - k_1 \cdot A \cdot Y - k_2 \cdot X \cdot Y + f \cdot k_5 \cdot A \cdot Z \\[2ex]&\frac{\mathrm{d}Z}{\mathrm{d}t} = k_3 \cdot A \cdot X - k_5 \cdot A \cdot Z.\end{align*}$

Therein, $X$ , $Y$ , and $Z$ are the concentrations of relevant intermediates of this reaction network, i.e., $\mathrm{HBrO}_2$ , $\mathrm{Br}^-$ , and Ce(IV). $A$ is the concentration of $\mathrm{BrO}_3^-$ , a constant. All $k_k$ and $f$ are parameters of the kinetics. For a more thorough discussion of the underlying physics, please consult the literature on the oregonator. We want to focus on the model implementation here.

Model workflow

Notation

We beginn by setting up the notation. Therefore, perform the following steps:

Go to the Notation tab in the “Model” section
Add a description of your notation
Add the base names below by clicking on “+ Name” at the bottom, entering these names, adding a description of the variable and confirming it
Go to the Indicies tab and add the index below
Once you have added all variables and symbols, click on “Save” at the top of the editor and save your notation at a suitable location in your folder structure

Base names

$A$ , concentration of $\mathrm{BrO}_3^-$ in mol/L
$f$ , Stoichiometric factor
$k$ , kinetic factor in various units
$t$ , time in s
$X$ , concentration of $\mathrm{HBrO}_2$ in mol/L
$Y$ , concentration of $\mathrm{Br}^-$ in mol/L
$Z$ , concentration of Ce(IV) in mol/L

Indices

$k$ , Index of kinetic factors 1..NK

The notation of this example has ID 182866.

Equations

Go to the Equation tab. Proceed as follows:

Load the notation you just created
Add an informative description of the current equation
Type the first model equation in LaTeX style into the window MosaicTex. For the differential operator, you must use \diff{Var}{IndependentVariable}
Save the equation at a suitable location in your folder structure
Click on “New” at the top of the editor and repeat the steps 1 and 2
Enter the second model equation
Save this second equation as well

The three equations are available in MOSAICmodeling with the IDs 182867, 182868, and 182869.

Equation system

Now, we will combine our equations to an equation system. To this end, go to the Equation System tab and carry on as described below:

Load your notation from above
Add a suitable description
Go to the tab Connected Elements and click on “+ EQU/EQS” at the bottom of the editor
In the opening popup window, select your three equations and confirm (you can also add them one-by-one). They will then appear in the overview of connected elements
Save your equation system

The equation system has ID 182870 within MOSAICmodeling.

Simulation workflow

Start by going to the “Simulation” section of MOSAICmodeling.

Equation system

Go to the tab Equation System and select your newly saved equation system from above.

Indexing

Proceed to the tab Indexing. Specify the maximum value of the index for the kinetic factors, i.e., 5. Click on Confirm Index Data.

Specification

Go to the tab Specifications and then select the tab Variables. Proceed as follows:

Add a helpful description of the simulation
In the specification window, you can see all possible variable classifications. Select ALL VARIABLES at the top. You should now see the four variables that are part of your model
Change the entry in the column Type for $A$ , $f$ and $k_1$ – $k_5$ to Design Value
Change the type classification of the variables $X$ , $Y$ , and $Z$ to state variable. The variable $t$ should already be classified as differential value. Your degree of freedom should now be zero
Assign the values given in Table 1 to the four variables
On the right side, set 2000.0 as end point for the differential variable
Save your variable specification by clicking on the “Save” button above the calculation of the degrees of freedom
Save your simulation by clicking on “Save” at the top of the editor

The simulation is available with ID 182871. The associated variable specification has ID 182872. Note that the values specificed for the state variables are used as initial conditions.

Name	Description	Value / Initial condition
$A$	Concentration of $\mathrm{BrO}_3^-$ in mol/L	$0.06$
$f$	Stoichiometric factor	$1.0$
$k_{k=1}$	Kinetic parameter	$1.34$
$k_{k=2}$	Kinetic parameter	$1.6 \cdot 10^9$
$k_{k=3}$	Kinetic parameter	$8.0 \cdot 10^3$
$k_{k=4}$	Kinetic parameter	$4.0 \cdot 10^7$
$k_{k=5}$	Kinetic parameter	$1.0$
$X$	Concentration of $\mathrm{HBrO}_2$ in mol/L	$7.8 \cdot 10^{-9}$
$Y$	Concentration of $\mathrm{Br}^-$ in mol/L	$8.502 \cdot 10^{-8}$
$Z$	Concentration of Ce(IV) in mol/L	$0.0190464$
$t$	Differential variable time in s	$0$

Table 1: Overview of design values and initial conditions.

Evaluation

Finally, we want to determine the solution of our model:

Go to the Evaluation tab, then go to the Code Generation tab
Select your preferred language specification, e.g., Matlab ODE std
Click on Generate Code
Go to the tab View Code and copy the code into your respective modeling environment, e.g., Matlab
Execute the code
You will notice that (i) the solution requires a significant amount of time, and (ii) the solution does not fluctuate at all. The long computational time is caused by ODE45, which is not suitable for stiff systems such as the ODE system in this example. Therefore, replace the following line with the one given below that

% [X,Y] = ode45(@(X,Y)calculateDifferentials(X,Y,PARAMS),X_INTERVAL,Y_INIT');
[X,Y] = ode23s(@(X,Y)calculateDifferentials(X,Y,PARAMS),X_INTERVAL,Y_INIT');

You will now receive a solution must faster, but without any fluctuation. The problem here is the accuracy of the applied solver. Now, change the respective options as shown below:

% [X,Y] = ode45(@(X,Y)calculateDifferentials(X,Y,PARAMS),X_INTERVAL,Y_INIT');
% [X,Y] = ode23s(@(X,Y)calculateDifferentials(X,Y,PARAMS),X_INTERVAL,Y_INIT');
[X,Y] = ode23s(@(X,Y)calculateDifferentials(X,Y,PARAMS),X_INTERVAL,Y_INIT',odeset('RelTol',1e-7,'AbsTol',1e-7));

Execute the code again. You should now observe a fluctuation of the three concentrations. This illustrates that – although MOSAICmodeling (or any other modeling tool) can help you with setting up the model equations – it does not absolve you from knowing about numerical issues and how to address them in the tool of your preference.