|Home | About | Journals | Submit | Contact Us | Français|
A method has been developed for parameter estimation using sensitivity analysis in a model for transport of substrates and ions between the blood stream and myocardial cells. Experimental data include indicator-dilution curves recorded from the coronary sinus of dogs and rabbits after bolus injections into the aorta. In addition, the intramyocardial density of deposition of radioactive microspheres was used to provide estimates of the heterogeneity of regional flows in the heart. The sensitivity analysis employed a calculation of the domains of dominance (regions of maximal impact defined with respect to the form of the indicator-dilution curves and the time after injection) for each parameter on the basis of their sensitivity functions.
Parameter estimation in a new two-barrier convection-diffusion model of coronary circulation has been achieved with great efficiency. Computation time was reduced to several minutes on an Interdata minicomputer as compared with 45 minutes to several hours using conventional optimization techniques such as STEPT (a direct search routine) or Simplex. The new parameter estimation method based on selective sensitivity is recommended for the class of problems which involve parameters with a well expressed distribution over the range of the independent variable.
Estimation of parameters of a complex nonlinear model by conventional objective function minimization methods (such as Marquardt, Fletcher-Powell, direct search or Simplex techniques) may be tedious. A feature common to all algorithms mentioned above is that they are designed primarily for minimization of an objective function, usually the sum of squares, giving optimal parameter estimates as by-products. Can a more direct process be generalized? When the influences of parameters are distributed over a range of independent variable, the sensitivity analysis provides information allowing a natural direct estimation procedure, and one which does not require the use of any overall objective function. That is especially useful in cases when the range of the dependent variable includes several orders of magnitude and the weighting functions are not known.
In general, sensitivity analysis concerns the effects of parameter variations on the behavior of dynamic systems [6, 7], for example, the stability and control properties of systems. Although it was strongly recommended  that before any optimization is attempted the sensitivity functions be plotted and examined for linear dependence in order to ensure the identifiability of parameters, only a few workers have used these for parameter estimation .
Consider a general functional relationship
where t is an independent variable, ŷ will be tabulated against t1, t2, …, tN; and = (p1, p2, …, pJ) a J-dimensional vector of parameters. The sensitivity function is
where Δŷj = ŷ (t, p1, …, pj + Δpj, …, pJ) − ŷ (t, p1, …, pj, …, pJ) and pj + Δpj is a small perturbation of the parameter pj. The impact function , a measure of the change in ŷ caused by perturbation of pj, is:
omitting the smaller terms of the Taylor series expansion. The domain of dominance of the parameter pj against the parameter ps is formally defined as the range of t where
Let ρ = Δpj/Δps. Then, for example, [− ρ, ρ] is the domain of dominance of the intercept p1 for a straight line ŷ =p1 + p2t. For an exponential ŷ=p1 + exp (p2t), the domain of dominance for p1 is given by [− ρ/p1 < t < ρ/p1]. For a Gaussian function , the domain of dominance of the mean p1 is given by the double range [p1 + w < t <p1 − v] and [p1 − w < t <p1 + v], where , the rest being the domain of dominance of the variance p2. These examples illustrate the existence of selective sensitivity: the dependent variable is more sensitive to perturbations of a parameter within its domain. Given the model function, ŷ, and the experimental curve, y, the difference ŷ − y at some maximal impact point tm for each parameter within its domain can be expanded into a linearized recurrence equation for iteration k:
The resulting system of linear algebraic equations can be solved by usual methods like Gaussian elimination to yield new parameter estimates. The use of a localized signed “distance” between model y and experimental y curves, as in Eq. (5), was described by Bassingthwaighte . In analyzing experimental functions subject to noise, Eq. (5) is useful only when the sensitivity function rises above the noise level. For Gaussian noise proportional to ŷ and with a standard deviation equal to β times the mean, parameter pj can be best estimated from (5) if the fractional error is greater than β/5:
at the maximal impact point tm, where 5 is an empirical choice.
The following algorithm was designed for parameter estimation using the sensitivity analysis: For each iteration k, (1) compute the matrix of J sensitivity functions (Eq. 2) at N time points; (2) for each parameter pj (j = 1, 2, …, J), find the domain of its dominance against all other parameters (Eq. 4); (3) for each parameter pj within its domain, find the maximal impact point tm at which the impact function (Eq. 3) has its maximum; (4) if the fractional error (Eq. 6) is less than β/5, remove the parameter pj from consideration for iteration k and go to step (2), but if all parameters are excluded, end;(5) calculate the new parameter estimate (Eq. 5); (6) if exceeds either a lower or an upper bound, set it equal to the bound; (7) if any of the termination criteria is satisfied, end; (8) set k = k + 1 and go to step (1).
The program terminates when either one of the following criteria is satisfied: (a) if (6) is violated; (b) if the number of iterations exceeds some preassigned maximum; or (c) if Sk < S, where S is some measure of fit. The algorithm requires as little as J + 1 model function evaluations at each time point per iteration. Domains can be manipulated by imposing some windows or gates along t-axis for screening the impact functions (3). A window for pj is defined as some limited range of t in which the impact function is computed while being set to zero outside this range. This makes it possible to maximize the rate of convergence for each parameter.
The model is a convection-diffusion model for a system composed of a set of similar parallel independent pathways having a heterogeneity of flows but identity of all other parameters. The capillary system impulse response, hc(t), for capillary-tissue model (see Fig. 1) can be represented as
where hi is the model, wi(fi) is the weighting function of the regional flow, fi, in each region i, t is time, and is the vector of governing parameters, as shown on Fig. 2. For a given solution, p includes (1) the ratio γ of interstitial, VISF, to capillary, Vc, volumes of distribution; (2) the capillary permeability-surface area product per unit flow, PSc/F; (3) the ratio θ of intracellular volume of distribution, VM, to Vc; and (4) the permeability-surface area product per unit flow for tissue cell membrane, PSM/F. All four parameters are dimensionless.
Assuming that the transit time through the large vessels is independent of the local capillary flows and that the system is linear and stationary, the large vessel transport function hLV(t) can be convoluted with hc(t) to produce the total system response, h(t):
where the asterisk denotes the convolution integral. hLV(t) is considered here to be the result of the convolution of the arterial and venous transport functions, i.e.
(see Fig. 1). For permeating tracers, the capillary-tissue exchange function hi of Eq. (7) is found as a solution to partial differential equations based on matter conservation. This solution may be numerical or in the form of an explicit function .
A series of experiments has been performed with different diffusible tracers (14C-D-glucose, 3H-2-deoxyglucose, and 14C-L-glucose) with a nonpermeating 131I human serum albumin serving as a vascular reference. The experimental techniques were similar to those described by Guller et al. . The primary objective of this study was to obtain reliable estimates of the essential parameters. During the preliminary stage of investigation it was noticed that the parameters of interest were not uniform with respect to their influence on the response: PSc and γ affect the upslope and the peak of the outflow curve, while PSM and θ tend to influence mainly the downslope and tail (tracer returning from cells), a separability of influences which was confirmed by sensitivity analysis and used in the optimization technique.
A set of test model curves with preassigned parameter values was produced and the parameters (γ, PScap, θ, PScell) were optimized from different starting points using three methods: the sensitivity algorithm (with and without windows), and STEPT (a direct search routine written by J.P. Chandler and distributed by the Quantum Chemistry Program Exchange, Oklahoma State University). The results are summarized in Table 1. It is seen that the windowed sensitivity estimation is about three times as fast as the no-window approach, and more than 10 times as fast as the direct search procedure (see Fig. 3).
The typical reference and diffusible test curves as well as the model sensitivity functions and their impact functions at the best fit are plotted in Fig. 4, which shows a clear distribution of parameter sensitivities over the time domain. At the next stage, a random noise was superimposed on the test curves using a pseudorandom number generator combined with a Gaussian (0, 1) distribution according to a scheme ŷ = y + εσ where ε was normally distributed with mean zero and variance one, and σ was set proportional to response, i.e. σ = βŷ for some 0 < β < 1. The results of parameter estimation for different values of β are given in Table 2. The maximal noise level most likely (with probability 0.95) does not exceed 1.96 β. A repeated application of noise to the test curve allowed us to obtain the confidence intervals for parameter estimates. 95% intervals at β = 0.05 for the parameter vector (γ, PSc, θ, PSM) are found to be γ ± 2.34, PSc ± 0.18, θ ± 3.5, PSM ± 0.1.
Finally, a set of experimental washout curves was subject to parameter estimation. The typical results are represented in Table 3. It was suggested that L-glucose does not enter the cell. Indeed, optimization with PSM fixed at zero produced a good fit (CV = 0.067) with the parameter values (γ = 8.74, PSc = 0.697) for the rabbit data set of Table 3. An illustration of the best fit to the experimental curve of column (3), Table 3, is given in Fig. 5.
It is our belief that sensitivity oriented parameter optimization can be used for a wide range of mathematical models. Generalizations to the multivariable multiresponse case are straightforward. Moreover, Eq. (1) can be perceived as a solution of some dynamic system of the form F(ÿ, , y, t, p) = 0, where dots denote derivatives with respect to time. In this case sensitivity functions can be computed from the linear differential equation (Tomovic, 1966):
The main advantage of sensitivity optimization is that it does not require casting specific objective functions (it is uncertain, for example, whether a weighted least squares, regular or logorithmic coefficient of variation or any other measure of fit is most suited for estimation of parameters in our model).
This work represents studies undertaken in partial fulfillment of the requirements of a Ph.D. program for M. Levin. The experiments were performed with the help of Mr. R. B. King. Hedi Nurk assisted with the illustrations. Supported by NIH grants HI 19139, HL 19135 and by Reynolds Industries.