# A template for parameter estimation with Matlab Optimization Toolbox; including dynamic systems

**1. Curve fitting**

A *weighted least squares fit* for a model which is
less complicated than the system that generated the data (a case of
so-called ‘undermodeling’).

System:

Data (noise model):

with additive zero mean white noise ZMWN with a standard
deviation of 0.05.

Model to fit the data:

with parameters to be fitted *q* = [*a,b*].

The model can be ‘forced’ to fit the second part of the data (1
< *x* < 2) by assigning different
weights to the data. Error function:

with *w*1 = 1 and *w*2
= 10.

To implement and solve the weighted least squares fitting problem in
Matlab the function LSQNONLIN of the Optimization Toolbox is used.
Optimization algorithms (in fact a minimization is performed) require
the user to specify an initial guess *q*0 for the
parameters. Here we use *q*0 = [0.1, -1 ].

Executing the Matlab program myFit.m results in:

`parameter_hat =`

`
3.1001 -1.9679`

However, here we do not obtain any information about the accuracy of the estimated . For example, we do not know how critical the fit depends on all digits in the parameter values returned by Matlab.

**2. Parameter estimation for a dynamic model**

In the second example we consider a dynamical system. If blood plasma
and a tissue or organ of interest can be considered as connected
compartments then the following model can be used to describe tissue
perfusion:

in which the arterial plasma concentration *Ca*(*t*)
is the input and the tissue concentration *Ct*(*t*)
is the output.

This linear differential equation model can be transformed into 2
descriptions often used in system and control engineering.

a) State-space format:

,

b) Transfer function:

The Control Toolbox from Matlab can be used to implement and
simulate this model.

The parameters *q* = [*Ktrans*, *u**e*]
need to be estimated from clinical data. Tissue perfusion can be
measured with Dynamic Contrast Enhanced MRI (DCE-MRI). A contrast agent
is injected into the bloodstream, circulates through the body and
diffuses into tissues and organs. The kinetics of contrast enhancement
is measured in a voxel in the tissue of interest, which can be
translated into a concentration. A dataset
is available.

The following noise (error) model is assumed/proposed:

with xi* *ZMWN and therefore the sum of
squared errors is used as estimator.

In this case the purpose of the model is not just to be able to
describe the data, but differences in the estimated parameter values
have a biological meaning and might be used for clinical diagnosis and
decisions. Hence it is critically important to assess the accuracy of
the estimates. For this problem the estimator is the Maximum Likelihood
Estimator (MLE) and it is possible to calculate the covariance matrix
of the parameter estimates, which quantifies the accuracy of the
estimate. (The inverse of the covariance matrix is known as the Fisher
Information Matrix.)

In the Matlab code compartment.m
it is shown how the covariance matrix can be calculated from the
outputs provided by the LSQNONLIN function. As initial guesses *q*0
= [2e-3, 0.1] will be used.

When the program is executed, first a figure is returned with
the input data (concentration in blood plasma, AIF, in blue), the
output data (tissue compartment, Ct, in red) and the model simulation
given the initial parameter values *q*0 (solid
black line).

Secondly the results of the parameter estimation (vector of estimated parameters, estimated parameter covariance matrix, and the parameter estimates with their relative standard deviation, in (%)):

`p =`

`
0.0097 0.2018`

`varp =`

` 1.0e-005 *`

`
0.0133 -0.0025
-0.0025 0.2574`

`Ktrans: 0.0096992 +/- 3.764%
ve: 0.20176 +/- 0.79515%`

Finally, a figure with a plot of the identified model and the data (top panel), and a plot of the difference between model and data (residuals, bottom).

**3. Template **

Finally a template
is provided to estimate a *subset* of the
parameters in a model (some parameters are assumed to be known and
therefore are fixed) and the model is composed of a set of coupled
first order *nonlinear differential equations*
(simulated with one of the Matlab ode-solvers).

An example in which this template is applied is a well-known compartment model for whole-body glucose homeostasis. Measured time-courses data of blood plasma insulin and glucose are used as input and output data, respectively, to estimate parameters that represent insulin sensitivity of a subject.

List of files:

- Van Riel, N.A.W. (2011) Template for parameter estimation with Matlab Optimization Toolbox; including dynamic systems DOI: 10.13140/2.1.3335.5846
- myFit.m
- compartment.m
- dataset
- template
- Van Riel, N.A.W. (2006) Parameter estimation in non-equidistantly sampled nonlinear state space models; a Matlab implementation