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

pdf

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:

                                                                            equation1                                                                 

Data (noise model):

                                                                              equation2                                                                   

with additive zero mean white noise ZMWN with a standard deviation of 0.05.
Model to fit the data:

                                                                                        equation3                                                                             

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:

                                        equation4                             

with w1 = 1 and w2 = 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 q0 for the parameters. Here we use q0 = [0.1, -1 ].

Executing the Matlab program myFit.m results in:

parameter_hat =

    3.1001   -1.9679

curvefit

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, ue] 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 q0 = [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 q0 (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: