# Using AMICI’s MATLAB interface¶

In the following we will give a detailed overview how to specify models in MATLAB and how to call the generated simulation files.

Note

The MATLAB interface requires the MathWorks Symbolic Math Toolbox for model import (but not for model simulation).

The Symbolic Math Toolbox requirement can be circumvented by performing model import using the Python interface. The resulting code can then be used from Matlab (see Compiling a Python-generated model).

Warning

Due to changes in the Symbolic Math Toolbox, the last MATLAB release with working AMICI model import is R2017b (see https://github.com/AMICI-dev/AMICI/issues/307).

## Specifying models in Matlab¶

This guide will guide the user on how to specify models in MATLAB.
For example implementations see the examples in the `matlab/examples`

directory.

### Header¶

The model definition needs to be defined as a function which returns a
`struct`

with all symbolic definitions and options.

```
function [model] = example_model_syms()
```

### Options¶

Set the options by specifying the respective field of the model struct

```
model.(fieldname) = value
```

The options specify default options for simulation, parametrisation and compilation. All of these options are optional.

field |
description |
default |
---|---|---|

.param |
default parametrisation ‘log’/’log10’/’lin’ |
‘lin’ |

.debug |
flag to compile with debug symbols |
false |

.forward |
flag to activate forward sensitivities |
true |

.adjoint |
flag to activate adjoint sensitivities |
true |

When set to `false`

, the fields `forward`

and `adjoint`

will speed up the
time required to compile the model but also disable the respective sensitivity
computation.

### States¶

Create the respective symbolic variables. The name of the symbolic variable can be chosen arbitrarily.

```
syms state1 state2 state3
```

Create the state vector containing all states:

```
model.sym.x = [ state1 state2 state3 ];
```

### Parameters¶

Create the respective symbolic variables. The name of the symbolic variable can
be chosen arbitrarily. Sensitivities will be derived for all *parameters*.

```
syms param1 param2 param3 param4 param5 param6
```

Create the parameters vector

```
model.sym.p = [ param1 param2 param3 param4 param5 param6 ];
```

### Constants¶

Create the respective symbolic variables. The name of the symbolic variable
can be chosen arbitrarily. Sensitivities with respect to *constants* will not
be derived.

```
syms const1 const2
```

Create the constants vector

```
model.sym.k = [ const1 const2 ];
```

### Differential equations¶

For time-dependent differential equations you can specify a symbolic variable
for time. This **needs** to be denoted by `t`

.

```
syms t
```

Specify the right hand side of the differential equation `f`

or `xdot`

```
model.sym.xdot(1) = [ const1 - param1*state1 ];
model.sym.xdot(2) = [ +param2*state1 + dirac(t-param3) - const2*state2 ];
model.sym.xdot(3) = [ param4*state2 ];
```

or

```
model.sym.f(1) = [ const1 - param1*state1 ];
model.sym.f(2) = [ +param2*state1 + dirac(t-param3) - const2*state2 ];
model.sym.f(3) = [ param4*state2 ];
```

The specification of `f`

or `xdot`

may depend on states, parameters and
constants.

For DAEs also specify the mass matrix.

```
model.sym.M = [1, 0, 0;...
0, 1, 0;...
0, 0, 0];
```

The specification of `M`

may depend on parameters and constants.

For ODEs the integrator will solve the equation \(\dot{x} = f\) and for DAEs the equations \(M \cdot \dot{x} = f\). AMICI will decide whether to use CVODES (for ODEs) or IDAS (for DAEs) based on whether the mass matrix is defined or not.

In the definition of the differential equation you can use certain symbolic
functions. For a full list of available functions see
`src/symbolic_functions.cpp`

.

Dirac functions can be used to cause a jump in the respective states at the specified time-point. This is typically used to model injections, or other external stimuli. Spline functions can be used to model time/state dependent response with unknown time/state dependence.

### Initial Conditions¶

Specify the initial conditions. These may depend on parameters on constants
and must have the same size as `x`

.

```
model.sym.x0 = [ param4, 0, 0 ];
```

### Observables¶

Specify the observables. These may depend on parameters and constants.

```
model.sym.y(1) = state1 + state2;
model.sym.y(2) = state3 - state2;
```

In the definition of the observable you can use certain symbolic functions.
For a full list of available functions see `src/symbolic_functions.cpp`

.
Dirac functions in observables will have no effect.

### Events¶

Specifying events is optional. Events are specified in terms of a trigger
function, a bolus function and an output function. The roots of the trigger
function defines the occurrences of the event. The bolus function defines the
change in the state on event occurrences. The output function defines the
expression which is evaluated and reported by the simulation routine on every
event occurrence. The user can create events by constructing a vector of
objects of the class `amievent`

.

```
model.sym.event(1) = amievent(state1 - state2,0,[]);
```

Events may depend on states, parameters and constants but *not* on observables.

For more details about event support see:

Fröhlich, F., Theis, F. J., Rädler, J. O., & Hasenauer, J. (2017). Parameter estimation for dynamical systems with discrete events and logical operations. Bioinformatics, 33(7), 1049-1056. doi:10.1093/bioinformatics/btw764.

### Standard deviation¶

Specifying standard deviations is optional. It only has an effect when computing adjoint sensitivities. It allows the user to specify standard deviations of experimental data for observables and events.

Standard deviation for observable data is denoted by `sigma_y`

```
model.sym.sigma_y(1) = param5;
```

Standard deviation for event data is denoted by `sigma_t`

```
model.sym.sigma_t(1) = param6;
```

Both `sigma_y`

and `sigma_t`

can either be a scalar or of the same dimension
as the observables / events function.
They can depend on time and parameters but must not depend on the states or
observables. The values provided in `sigma_y`

and `sigma_t`

will only be used
if the value in `D.Sigma_Y`

or `D.Sigma_T`

in the user-provided data struct is
`NaN`

. See simulation for details.

### Objective Function¶

By default, AMICI assumes a normal noise model and uses the corresponding negative log-likelihood

as objective function. A user provided objective function can be specified in

```
model.sym.Jy
```

As reference see the default specification of `this.sym.Jy`

in `amimodel.makeSyms`

.

## SBML¶

AMICI can also import SBML models using the command `SBML2AMICI`

.
This will generate a model specification as described above, which may be
edited by the user to apply further changes.

## Model Compilation¶

The model can then be compiled by calling `amiwrap.m`

:

```
amiwrap(modelname,'example_model_syms',dir,o2flag)
```

Here `modelname`

should be a string defining the name of the model, `dir`

should be a string containing the path to the directory in which simulation
files should be placed and `o2flag`

is a flag indicating whether second order
sensitivities should also be compiled.
The user should make sure that the previously defined function
`example_model_syms`

is in the user path. Alternatively, the user can also
call the function `example_model_syms`

```
[model] = example_model_syms()
```

and subsequently provide the generated struct to `amiwrap(...)`

, instead of
providing the symbolic function:

```
amiwrap(modelname,model,dir,o2flag)
```

In a similar fashion, the user could also generate multiple models and pass
them directly to `amiwrap(...)`

without generating respective model
definition scripts.

### Compiling a Python-generated model¶

For better performance or to avoid the Symbolic Math Toolbox requirement,
it might be desirable to import a model in Python and compile the
resulting code into a mex file. For Python model import, consult the
respective section of the Python documentation. Once the imported
succeeded, there will be a `compileMexFile.m`

script inside the newly
created model directory which can be invoked to compile the mex file.
This mex file and `simulate_*.m`

can be used as if fully created by
matlab.

#### Using Python-AMICI model import from Matlab¶

With recent matlab versions it is possible to use the AMICI python package from within Matlab. This not quite comfortable yet, but it is possible.

Here for proof of concept:

Install the python package as described in the documentation

Ensure

`pyversion`

shows the correct python version (3.6 or 3.7)Then, from within the AMICI

`matlab/`

directory:

```
sbml_importer = py.amici.SbmlImporter('../python/examples/example_steadystate/model_steadystate_scaled.xml')
sbml_importer.sbml2amici('steadystate', 'steadystate_example_from_python')
model = py.steadystate.getModel()
solver = model.getSolver()
model.setTimepoints(linspace(0, 50, 51))
rdata = py.amici.runAmiciSimulation(model, solver)
result = struct(py.dict(rdata.items()))
t = double(py.array.array('d', result.ts))
x = double(py.array.array('d', result.x.flatten()))
x = reshape(x, flip(double(py.array.array('d', result.x.shape))))
plot(t, x)
```

## Model simulation¶

After the call to `amiwrap(...)`

two files will be placed in the specified
directory. One is a `_modelname_.mex`

and the other is `simulate_*modelname*.m`

.
The mex file should never be called directly. Instead the MATLAB script, which
acts as a wrapper around the .mex simulation file should be used.

The `simulate_ _modelname_.m`

itself carries extensive documentation on how to
call the function, what it returns and what additional options can be
specified. In the following we will give a short overview of possible function
calls.

### Integration¶

Define a time vector:

```
t = linspace(0,10,100)
```

Generate a parameter vector:

```
theta = ones(6,1);
```

Generate a constants vector:

```
kappa = ones(2,1);
```

Integrate:

```
sol = simulate_modelname(t,theta,kappa,[],options)
```

The integration status will be indicated by the `sol.status`

flag. Negative
values indicated failed integration. The states will then be available as `sol.x`

.
The observables will then be available as `sol.y`

. The event outputs will then
be available as `sol.z`

. If no event occurred there will be an event at the end
of the considered interval with the final value of the root function is stored
in `sol.rz`

.

Alternatively the integration can also be called via

```
[status,t,x,y] = simulate_modelname(t,theta,kappa,[],options)
```

The integration status will be indicated by the flag `status`

. Negative
values indicated failed integration. The states will then be available as `x`

.
The observables will then be available as `y`

. No event output will be given.

### Forward Sensitivities¶

Set the sensitivity computation to forward sensitivities and integrate:

```
options.sensi = 1;
options.sensi_meth = 'forward';
sol = simulate_modelname(t,theta,kappa,[],options)
```

The integration status will be indicated by the `sol.status`

flag. Negative
values indicate failed integration. The states will be available as `sol.x`

,
with the derivative with respect to the parameters in `sol.sx`

.
The observables will be available as `sol.y`

, with the derivative with respect
to the parameters in `sol.sy`

. The event outputs will be available as `sol.z`

,
with the derivative with respect to the parameters in `sol.sz`

. If no event
occured there will be an event at the end of the considered interval with the
final value of the root function stored in `sol.rz`

, with the derivative with
respect to the parameters in `sol.srz`

.

Alternatively the integration can also be called via

```
[status,t,x,y,sx,sy] = simulate_modelname(t,theta,kappa,[],options)
```

The integration status will be indicated by the status flag. Negative values
indicate failed integration. The states will then be available as `x`

, with
derivative with respect to the parameters in `sx`

. The observables will then
be available as `y`

, with derivative with respect to the parameters in `sy`

.
No event output will be given.

### Adjoint sensitivities¶

Set the sensitivity computation to adjoint sensitivities:

```
options.sensi = 1;
options.sensi_meth = 'adjoint';
```

Define Experimental Data:

```
D.Y = [NaN(1,2)],ones(length(t)-1,2)];
D.Sigma_Y = [0.1*ones(length(t)-1,2),NaN(1,2)];
D.T = ones(1,1);
D.Sigma_T = NaN;
```

The `NaN`

values in `Sigma_Y`

and `Sigma_T`

will be replaced by the
specification in `model.sym.sigma_y`

and `model.sym.sigma_t`

. Data points
with `NaN`

value will be completely ignored.

Integrate:

```
sol = simulate_modelname(t,theta,kappa,D,options)
```

The integration status will be indicated by the sol.status flag. Negative
values indicate failed integration. The log-likelihood will then be available
as `sol.llh`

and the derivative with respect to the parameters in
`sol.sllh`

. Note that for adjoint sensitivities no state, observable and
event sensitivities will be available. Yet this approach can be expected to be
significantly faster for systems with a large number of parameters.

### Steady-state sensitivities¶

This will compute state sensitivities according to the formula

In the current implementation this formulation does not allow for conservation laws as this would result in a singular Jacobian.

Set the final timepoint as infinity, this will indicate the solver to compute the steadystate:

```
t = Inf;
```

Set the sensitivity computation to steady state sensitivities:

```
options.sensi = 1;
```

Integrate:

```
sol = simulate_modelname(t,theta,kappa,D,options)
```

The states will be available as `sol.x`

, with the derivative with respect
to the parameters in `sol.sx`

. The observables will be available as `sol.y`

,
with the derivative with respect to the parameters in `sol.sy`

. Notice that
for steady state sensitivities no event sensitivities will be available. For
the accuracy of the computed derivatives it is essential that the system is
sufficiently close to a steady state. This can be checked by examining the
right hand side of the system at the final time-point via `sol.diagnosis.xdot`

.