PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Med Prog Technol. Author manuscript; available in PMC 2010 August 30.
Published in final edited form as:
Med Prog Technol. 1980 June; 7(2-3): 119–124.
PMCID: PMC2929982
NIHMSID: NIHMS204555

Sensitivity Analysis in Optimization of Time-Distributed Parameters for a Coronary Circulation Model

Abstract

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.

Keywords: Parameter estimation, Sensitivity analysis, Myocardial transport, Heterogeneity of flow, Transcoronary circulation

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.

Theory

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 [2] 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 [4].

Consider a general functional relationship

y^=y^(t,p)
(1)

where t is an independent variable, ŷ will be tabulated against t1, t2, …, tN; and [p with tilde] = (p1, p2, …, pJ) a J-dimensional vector of parameters. The sensitivity function is

y^Pjy^/pj=Δy^j/Δpj
(2)

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 Δy^j2, a measure of the change in ŷ caused by perturbation of pj, is:

Δy^j2=y^Pj2Δpj2
(3)

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

Δy^j2>Δy^s2
(4)

Let ρ = Δpjps. 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 y^=exp[12(tp1P2)]/p22π, the domain of dominance of the mean p1 is given by the double range [p1 + w < t <p1v] and [p1w < t <p1 + v], where 0>w=p2(ρρ2+4)/2,v=p2(ρ+ρ2+4)/2>0, 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:

y^(tm,p+Δp(k))y(tm)=j=1Jy^Pj(pj(k)=pj(k+1))
(5)

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 [1]. 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:

|y^(tm,p(k))y(tm)y^(tm,p(k))|>β/5
(6)

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 pj(k+1) 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

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

Fig. 1
The capillary-time model. For explanation of symbols, see text

hc(t)=i=1nwifihi(t,fi,p)
(7)

where hi is the model, wi(fi) is the weighting function of the regional flow, fi, in each region i, t is time, and [p with tilde] 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.

Fig. 2
Model for mass transport through an organ with heterogeneous regional flow. The asterisks denote convolution

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):

h(t)=hLV(t)hc(t)
(8)

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.

hLV(t)=hA(t)hv(t)
(9)

(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 [5].

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. [3]. 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.

Results

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).

Fig. 3
The comparative convergence pattern for STEPT (a) sensitivity optimization without (b) and with windows (c). N is the number of model function computations; R is current parameter value as a fraction of final value
Table 1
Parameter optimization for test curves with known parameter values using STEPT and sensitivity optimizing routines. Parameter vector consisted of γ, PScap, θ, and PScell Elapsed time reported here is the total time spent in the system ...

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 α=maxt|y^yy^| 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.

Fig. 4
The typical test reference hR(t) and diffusible hD(t) curves (a), the model sensitivity (b) and impact functions (c) at the best fit for each parameter; on semilogarithmic scale
Table 2
Test curve no. 1 with different noise levels. Final values of parameter estimates are averaged over 10 simulation runs

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.

Fig. 5
The typical best fit to experimental diffusible curve. Dog data, diffusible tracer: 3H-2-deoxyglucose. Coefficient of variation = 0.0314. Insert: distribution of relative regional flows for this data set. Fi = regional flow, [F with macron] = mean ...
Table 3
Typical parameter estimates from experimental data. Assumed maximal noise level is 10%

Discussion

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, 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):

Fy¨y.Pj+Fy.y.Pj+FyyPj+FPj=0
(10)

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).

Acknowledgments

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.

References

1. Bassingthwaighte JB. Plasma indicator dispersion in arteries of the human leg. Circ Res. 1966;19:332–246. [PMC free article] [PubMed]
2. Beck TJ, Arnold KT. Parameter estimation in engineering and science. Wiley; New York: 1977.
3. Guller B, Yipintsoi T, Orvis AL, Bassingthwaighte JB. Myocardial sodium extraction at varied coronary flows in the dog: Estimations of capillary permeability by residue and outflow detection. Circ Res. 1975;37:359–378. [PMC free article] [PubMed]
4. Howland TL, Vaillancourt R. A generalized curve-fitting procedure. SIAM. 1961;9:165–168.
5. Rose CP, Goresky CA, Bach GG. The capillary and sarcolemmal barriers in the heart. Circ Res. 1977;41:515–533. [PubMed]
6. Tomovic R. The Role of sensitivity analysis in engineering problems. Sensitivity methods in control theory. In: Radanovic L, editor. Proc IF AC Internat Symposium at Dubrovnik. Pergamon Press; New York: 1966.
7. Tomovic R, Vukobratovic M. General sensitivity theory. American Elsevier; New York: 1972.