where K1 is the rate constant of flow from vasculature to tissue, k2 is the rat constant of flow from tissue to vasculature, A(t) is the activity in the brain, and Ca(t) is the activity in the blood.
The MATLAB leastsq function provided with the Optimization Toolbox is the obvious choice for solving this problem, since it performs a non-linear least squares fit to a multi-dimensional problem. Please see the MATLAB documentation for a detailed description of how the leastsq function works, and what options it accepts. For our purposes, its default tolerances and fitting algorithm (Levenberg-Marquardt) will be used.
There are several steps to performing a curve fit in MATLAB:
function sim = simrat (vars, times, Ca) expfun = exp(-vars(2)*times); conv = nconv(Ca, expfun, times(2)-times(1)); conv = conv (1:length(times)); sim = vars(1)*conv;This short function returns the values of A(t) given the fitting variables (in this case, K1=vars(1), and k2=vars(2)), a vector of times, and a vector of values for Ca. The first line calculates the exponential part of the equation, the second line calculates the convolution, the third line clips the convolution to the length that we are interested in, and the fourth line sets the return value (K1 times the convolution).
function error = fitrat (vars, ts_even, plasma_even, brain_even) sim = simrat (vars, ts_even, plasma_even); error = brain_even - sim;
[final, options, f, j] = leastsq ('fitrat', [K1 k2], , , ... ts_even, plasma_even, brain_even);leastsq will call the function fitrat to get the residuals, and pass it [K1 k2] (which will be received in vars), ts_even, plasma_even, and brain_even in that order. The leastsq function passes back the final parameters (the solution) in final, the options vector (giving such information as the number of iterations), and both the residuals and the Jacobians at the solution.
[conf, var] = confint (final, f, j);This function takes the solution to the least squares problem (final), the residuals at the solution (f), and the Jacobians at the solution (j). It returns the 95% confidence interval and the variance (please see the MATLAB Optimization Toolbox documentation for further details, or type help confint in MATLAB).