- Email: [email protected]
winstodec : a stochastic deconvolution interactive program for physiological and pharmacokinetic systems
Computer Methods and Programs in Biomedicine 67 (2001) 67 – 77 www.elsevier.com/locate/cmpb
WINSTODEC:
a stochastic deconvolution interactive program for physiological and pharmacokinetic systems
Giovanni Sparacino a, Gianluigi Pillonetto a, Massimo Capello a, Giuseppe De Nicolao b, Claudio Cobelli a,* a
Dipartimento di Elettronica ed Informatica, Uni6ersita` di Pado6a, Via Gradenigo 6 /A, 35131 Pado6a, Italy b Dipartimento di Informatica e Sistemistica, Uni6ersita` di Pa6ia, Via Ferrata 1, 27100 Pa6ia, Italy Received 6 June 2000; received in revised form 15 November 2000; accepted 24 November 2000
Abstract Deconvolution allows the reconstruction of non-accessible inputs (e.g. hormone secretion rate) from their causally-related measurable effects (e.g. hormone plasma concentration). Deconvolution is challenging under several aspects both general (e.g. determination of a suitable trade-off between data fit and solution smoothness in order to contrast ill-conditioning, assessment of the confidence intervals) as well as specific of physiological systems (e.g. non-uniform and infrequent data sampling). Recently, a stochastic regularization approach has been proposed and validated to handle these difficulties (De Nicolao et al., Automatica 33 (1997) 851– 870). In this paper, an interactive program, WINSTODEC, is presented to allow the clinical investigator to easily obtain the solution of a deconvolution problem by this approach. © 2002 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Software; Inverse problems; Smoothing; Biomedical systems; Endocrinology; Regularization
1. Introduction In the quantitative study of physiological and pharmacokinetic systems, deconvolution allows the reconstruction of important but non-accessi
A preliminary version of this work was presented at the Fourth IFAC Symposium on Modeling and Control of Biomedical Systems, Karlsburg, Germany, 30 March–1 April 2000. * Corresponding author. Present address: Department of Electronics and Informatics, University of Padova, Via Gradenigo 6/A, 35100 Padova, Italy. Tel.: + 39-49-8277616; Fax: + 39-49-8277699. E-mail address: [email protected] (C. Cobelli).
ble inputs (e.g. the secretion rate of a gland, the rate of appearance of a drug into the plasma after an oral administration) from samples of their causally-related measurable effects (e.g. the time course of hormone or drug in the plasma), assuming the kinetics to be linear and the impulse response to be available (e.g. from population studies or suitably designed input –output experiments). Deconvolution is well known to be challenging under several aspects both general (e.g. determination of a suitable trade-off between data fit and solution smoothness in order to contrast ill-conditioning, estimation of the confidence intervals of the deconvoluted input) as well as spe-
0169-2607/02/$ - see front matter © 2002 Elsevier Science Ireland Ltd. All rights reserved. PII: S 0 1 6 9 - 2 6 0 7 ( 0 0 ) 0 0 1 5 1 - 6
68
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
cific of physiological systems (e.g. data sampling is often non-uniform and/or infrequent, the deconvoluted input usually must be non-negative, the kinetics can be time varying). Recently, a non-parametric stochastic deconvolution method has been proposed to handle the presented difficulties [1]. In particular, thanks to the concept of virtual grid, an input profile continuous to any desired order of time derivatives can be obtained even when the available data are collected on a sparse (possibly non-uniform) time schedule. Moreover, thanks to the Bayesian framework, regularization criteria based on the Maximum Likelihood principle and confidence intervals accounting for the bias error can be derived. Non-negativity constraints can also be included. Suitable numerical algorithms permit the quick and efficient solution of deconvolution problems in a wide variety of situations, including when the system is time varying. These algorithms are based on a singular value decomposition (SVD) strategy that, by putting the deconvolution problem in diagonal form, speeds up the trial and error procedure for the tuning of the regularization parameter that controls the balance between data fit and solution smoothness. By using a prototype software, developed within the Matlab environment, the performance of this algorithm was compared in Ref. [2] with that of other three algorithms/programs from the literature [3– 5] on about 50 benchmark problems, both real and simulated, including hormone secretion problems in spontaneous as well as stimulated conditions, drug rate of absorption problems and the ‘Hunt benchmark’ [6]. This wide platform of benchmark problems covered a wide spectrum of situations: many/a few data; frequent/infrequent sampling; high/mild ill-conditioning; bad/good signal-to-noise ratio; time variance/invariance of the dynamic system. The main outcome of the comparison in Ref. [2] was that, on average, the stochastic deconvolution approach of Ref. [1] proved to be, thanks to both its firm statistical basis and numerical efficiency of the SVD strategy, more robust than the other three algorithms considered. The aim of the present paper is to introduce the reader to the interactive program WINSTODEC
(where STODEC stands for STOchastic DEConvolution), the computational core of which is basically the same as that of the prototype employed in Ref. [2]. The major merit of WINSTODEC is the new graphical user interface, which should render it extremely easy to use, despite the sophistication of theory and algorithms, for both physiologists, physicians, and pharmacokineticists as well as biomedical engineers without specific expertise in deconvolution.
2. Methods This section provides the reader with a brief description of the theory behind the WINSTODEC software. Such an overview is not meant to be complete. A full and detailed documentation of the theory is reported in Ref. [1].
2.1. The ingredients of a decon6olution problem In a linear single-input single-output system, the input u(t) and the output c(t) are linked by the convolution integral: c(t)=
&
t
g(t,~)u(~) d~
(1)
−
where g(t,~) is the kernel of the system. If the system is time invariant, then g(t,~)= g(t− ~) and g(t) is the (unit) impulse response of the system. Independently of the deconvolution approach under adoption, the two main pieces of information required to estimate the system input via deconvolution are the samples of the output and the system kernel. For instance, consider the problem of reconstructing the secretion rate of C-peptide (equal to the secretion rate of insulin, equimolarly secreted from the pancreas) in a normal man in spontaneous conditions. The ingredients are depicted in Fig. 1. The top panel displays the (noisy) samples of the C-peptide concentration in plasma collected for nearly 20 h every 10 min (therefore, the output time-series consists of n=115 samples). The bottom panel depicts the sum of two exponentials that describes the impulse response g(t) in the same subject (solid line). The impulse response
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
69
was obtained by fitting the C-peptide plasma concentration samples (diamonds) following an intravenous pulse of biosynthetic C-peptide administered in the same subject. Usually, information on the precision of the plasma concentration data is also available, e.g. measurement errors have a constant variance (known or unknown) or a constant coefficient of variation (CV) (known or unknown). In this representative example, for instance, the CV of the C-peptide plasma concentration data is 5%.
extremely sensitive to noise in the data because of ill-conditioning. Then, the idea is to establish a trade-off between fidelity to the data and smoothness of the solution, measured by a suitable penalty functional. The input estimate is represented by an N-dimensional vector uˆ, the entries of which are the samples of the input over the virtual grid chosen by the user. If y denotes the observed vector (size n)
2.2. The stochastic regularization approach
where v is the (unknown) measurement error vector affecting the n samples of c(t) in Eq. (1), and G is the n× N ‘convolution matrix’ obtained by suitably discretizing the system impulse response. Then, the estimate uˆ can be computed as
WINSTODEC implements a regularization approach [7,8]. In fact, the most direct solution of the input estimation problem (i.e. the input that, once convoluted with the impulse response, guarantees perfect adherence to the output data) is
Fig. 1. The ingredients of a representative deconvolution problem (reconstruction of C-peptide secretion rate in a normal man). Top, C-peptide plasma concentration; bottom, C-peptide impulse response (solid curve) obtained by fitting a sum of two exponentials against C-peptide samples (diamonds) following an intravenous pulse of biosynthetic C-peptide administered in the same man (for both panels, time units are in minutes).
y = Gu +v
uˆ = (GTB − 1G+ k FTF) − 1GTB − 1y
(2)
(3)
where B is an n×n matrix equal, apart from a scale factor | 2, to the covariance matrix of v, F is an N× N penalty matrix, and k\ 0 is the socalled regularization parameter. In keeping with the theory of regularization, devised in the early 1960s by Phillips [7] and Tikhonov [8], matrix F is such that Fu is the vector of the mth time differences of u, m being a positive integer parameter. In other words, F= Z m, where Z is the square n-dimension lower-triangular Toeplitz matrix, the first column of which is [1, − 1, 0, …, 0]T. In the literature, there are only a few contributions aimed at theoretically addressing the choice of m (for example, Ref. [9]). A possible reason is that the deconvoluted profile is usually not very sensitive to the value of m. As a consequence, the typical choices for m are restricted to m= 1 or m= 2, depending on the particular problem (the user finds out a posteriori which is the best value between the two). As far as the regularization parameter k is concerned, its choice is a classical topic of the deconvolution literature, where many regularization criteria have been proposed in order to avoid subjectivity in its determination. In fact, k is a crucial parameter that governs the trade-off between data fit (measured by the squared norm of the residuals vector) and solution smoothness (measured by the energy of the mth-order time
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
70
derivatives of the input). The most known regularization criteria are heuristic or deterministic, e.g. Discrepancy/Twomey [10], Generalized CrossValidation [11], Minimum Risk [12]. The stochastic approach [1] allows the derivation of Maximum Likelihood (ML) criteria, usable when the variance (or CV) of the measurement error is either known or unknown. For instance, when the measurement error variance is known, k is iteratively adjusted until uˆ TFTFuˆ =
q(k) k
(4)
where q(k) =trace[B − 1G(GTB − 1G +kFTF) − 1GT]
(5)
are the so-called equivalent degrees of freedom. The stochastic approach also allows the estimation of the confidence intervals of the deconvoluted input from the diagonal of the covariance matrix of the estimation error e= u −uˆ : cov(e) =| 2(GTB − 1G+ kFTF) − 1
(6)
This feature represents a key advantage over the other deterministic regularization approaches previously proposed in the literature; see, for example, Ref. [6] where the difficulty of handling the bias error is explicitly admitted. The treatment of the non-uniform/infrequent sampling case relies on the introduction of the virtual grid. The virtual grid is a uniform N-size grid, chosen by the user, over which the input function is reconstructed. This grid is independent of the sampling grid (the only constraint being that the virtual grid must contain the sampling grid) and can be arbitrarily fine (i.e. N n). As the virtual grid becomes finer (i.e. N ), an input profile that is asymptotically continuous up to any desired order of time derivatives is eventually obtained, even when the available data are collected on a sparse (possibly non-uniform) time schedule. The degree of continuity reflects the value of the integer parameter m selected by the user in the definition of the penalty matrix F. For instance, the typical choices of m =2 or m = 1 lead to input functions asymptotically continuous at least up to the third or first time derivatives, respectively [1].
Suitable numerical algorithms permit the quick and efficient solution of deconvolution problems in a wide variety of situations. These algorithms are based on a SVD strategy that, diagonalizing the problem, speeds up the trial and error procedure for the determination of the regularization parameter. This will enable the WINSTODEC user to quickly compute solutions for several different values of k (corresponding to different regularization criteria). Finally, the regularization method implemented in WINSTODEC can be constrained for non-negativity. In this case, the solution does not exist in closed form and an iterative algorithm, the constrained conjugate gradient method, is employed to solve the problem: min(y−Gu)TB − 1(y−Gu) +kuTFTFu u]0
(7)
3. The software This section provides the user with an overall description of the software, the use of which will be exemplified in detail in Section 4. When WINSTODEC is launched, a Welcome Window opens that asks the user to supply the ingredients of the problem (i.e. data, system impulse response, …) and the regularization settings (i.e. regularization criterion, order of the penalty functional, …). This can be achieved by retrieving a file, previously stored on the disk, or interactively by filling in a sequence of windows. The interactive modality is that of interest in this paper. For instance, in the Output Window, the user can enter the data and the measurement error statistical description; in the System Window, the user can enter the impulse response in either a parametric or non-parametric format (i.e. amplitudes and exponents of a sum of decaying exponentials or a sequence of impulse response samples); in the Input Window, the user can specify name and units of the unknown input. Once this information has been supplied, the Regularization Tools Window opens, by which the user can specify the virtual grid over which the input must be reconstructed and the order of the penalty functional to be adopted (F in Eq. (3)).
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
Once this window has been closed, deconvolution can be invoked and results are placed in the Deconvolution Window. This is a key window in the program, since it enables the user to modify the settings of the algorithm and perform deconvolution again. In particular, the user can select from a menu different regularization criteria (Maximum Likelihood, Discrepancy/Twomey and Minimum Risk, when the measurement error variance is known, and Maximum Likelihood and Generalized Cross-Validation, when the measurement error variance is unknown) or force the regularization parameter to equal a value decided by himself/herself on a subjective basis (independent of any criterion). In addition, by operating on suitable switches, the Deconvolution Window allows the user to require the introduction of the non-negativity constraint and/or the calculation of the confidence intervals of the estimated input. The outcome of the deconvolution process consists of the estimated input profile, the reconvoluted profile and the residuals, the value of the regularization parameter together with the corresponding equivalent degrees of freedom and, when requested, the confidence intervals of the estimate. All these quantities are represented within the Deconvolution Window. Graphs of the signals are displayed by default, but the user has the possibility to invoke a window where tabular values can be inspected. A typical problem with n= 115 samples and N =1150 unknown input levels like the example of this paper is solved in less than 1 min on an ordinary 450 MHz CPU/64 Mb RAM PC (an additional minute is needed to estimate the confidence intervals). From the Deconvolution Window, the user has the possibility to access to a powerful Print Menu (which enables him/her to customize panels and labels of the plots) and save the numerical results on the disk together with the settings employed. If the deconvoluted profile is not satisfactory, from the Deconvolution Window, the user can change the settings and perform deconvolution again. For instance, he/she can change the regularization criterion, modify the penalty functional, require the confidence intervals, constrain the estimate to be non-negative. When deconvolution is repeated, the computation of the solution is usually very
71
Fig. 2. The Welcome Window allows the user to enter the ‘ingredients’ of a deconvolution problem (e.g. output samples, impulse response model) either in an interactive manner or by retrieving a file previously stored on the disk.
fast (a few seconds) since, in most cases, the changes (e.g. the regularization criterion of adoption) do not affect the previously computed SVD, which represents the most computationally-demanding step of the numerical algorithm. From the Deconvolution Window, the user can also invoke menus that allow him/her to graphically compare up to ten different solutions of the problem obtained at different times (e.g. by employing different regularization settings or impulse response models or measurement error statistics description). Once closed, a problem can be later retrieved from the Welcome Window, provided that it was saved on the disk. This feature allows the user to refine in any moment the solution of a problem (e.g. by asking the program to compute the confidence intervals or by using an improved impulse response model).
4. Example In this section, the use of WINSTODEC is exemplified by solving in steps the problem of Fig. 1. Once the program is launched, the Welcome Window opens (Fig. 2), by which a new problem can be defined. Here, we consider the case in which the data are supplied interactively. In the Welcome Window of Fig. 2, the fact that any of the
72
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
labels ‘Input’, ‘System’ and ‘Output’ appears in Italic fonts reminds the user that information is still to be supplied. By double-clicking on the Output arrow, the Output Window opens (Fig. 3), and the time-series of Fig. 1 (top) can be pasted in tabular format (time, value). Labels for name (‘C-peptide concentration’) and units (‘pmol/ml’ for the y axis and ‘minutes’ for the x axis) can be supplied. The ‘Plot’ button in the Output Window (see Fig. 3) allows the user to have an x – y scatter graph of the data entered (not shown). If a threecolumn table were supplied, the third column would be interpreted as the column of the measurement error standard deviations. When the third column is not specified, as in this example, the information on the measurement error statistics (here, CV constant equal to 5%) can be entered by a suitable menu that is automatically opened by the program. This menu (not shown) proposes five choices: (a) standard deviation constant with known value to be specified; (b) stan-
Fig. 3. The Output Window allows the user to enter samples, units and uncertainty of the measured output (the representative problem of Fig. 1 is hereafter considered).
Fig. 4. The System Window allows the user to enter either a parametric (sum of exponentials model) or non-parametric description of the impulse response.
dard deviation constant but unknown; (c) standard deviation dependent on the output through a generic law (e.g. quadratic) to be specified by the user; (d) coefficient of variation constant with known value to be specified by the user; (e) coefficient of variation constant but unknown. In this example, option (d) is selected. Once the Output Window has been filled in and closed (the font of the word ‘Output’ in the Welcome Window switches from the original Italics (see Fig. 2) to Roman), one can open the Input Window by double-clicking on the input arrow of Fig. 2. By this window (not shown), name and units of the unknown input (‘C-peptide secretion rate’, ‘pmol/min’) are entered. Once the Input Window is closed, the System Window can be opened (Fig. 4). In this window, the parameters of a sum of decaying exponentials impulse response model can be supplied interactively (in this example, amplitudes 8.947×10 − 5 and 2.999× 10 − 4, exponents 0.0228 and 0.197 for the C-peptide system of Fig. 1) or loaded from a file (e.g. one can have on the disk a library of population
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
impulse response models for a number of hormones). Alternatively, one could employ a nonparametric description of the impulse response (see the non-parametric description radio button in Fig. 4) and paste, in the window that would be opened, a two-column table containing the samples of the impulse response (not shown). These impulse response samples (which can also be loaded from a file) would be linearly interpolated in order to have a continuous-time description of g(t). In both cases, the ‘Plot’ button of the System Window (see Fig. 4) allows the user to graph the impulse response just entered (not shown). Having closed the System Window, the ‘Next’ command (see Fig. 2) becomes active. One can still revise any of the ingredients of the problem by double-clicking again on its symbol in the Welcome Window, but if ‘Next’ is double-clicked,
73
the Regularization Settings Window appears (not shown), where the virtual grid period and the order of the penalty functional can be chosen. The default values proposed by the software for the virtual grid period and the penalty functional order are, respectively, the greatest common divisor of the elements of the sampling grid and 2. In this particular C-peptide secretion problem, we let the virtual grid period be equal to 1 min (this allows the estimation of a ‘quasi-continuous’ input profile since it is approximated by a piecewise constant function with 1150 levels) and m be equal to 2 (a value that was seen to work slightly better than m= 1 in this particular problem). When the Regularization Settings Window is closed, the Deconvolution Window appears, where, after having selected from the proper menu a regularization criterion (here ‘Discrepancy’), de-
Fig. 5. The Deconvolution Window as it appears after having adopted the Discrepancy regularization criterion for the problem of Fig. 1.
74
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
Fig. 6. The Deconvolution Window as it appears after having adopted the Maximum Likelihood criterion to refine the solution of Fig. 5.
convolution can be launched. The results then appear in the Deconvolution Window in a graphical fashion (Fig. 5). In particular, three subplots are presented. The top subplot displays the deconvoluted input. The middle subplot depicts the fit, i.e. the reconvolution versus the data. The bottom subplot shows the residuals standardized by the standard deviation of the measurement errors. Interestingly, by suitably clicking on any of the subplots, one can inspect a zoom view, change the labels, the curve colour, and so on. At the right of the top subplot, the program reminds the user of the settings (i.e. m value and regularization criterion) and also reports the values achieved for k, q(k) and the weighted residuals sum of squares (WRSS) (here k= 6.6996 ×10 − 3, q(k) = 35.2466, WRSS =115.007). From the Deconvolution Window, a Print menu can be invoked (not shown),
from which one can decide how to arrange the subplots before printing the results on the printer or capturing the image. For instance, one can decide to visualize only the deconvolution profile and the reconvolution versus data (not shown). From the Deconvolution Window (Fig. 5), a ‘Save’ button is also available that allows the user to store the results on the disk with an arbitrary name (e.g. ‘solutionc 1.mat’). From the Deconvolution Window, the regularization settings can be modified and deconvolution can be launched again. For example, after having changed the regularization criterion of adoption to ‘Maximum Likelihood’ and clicked on the ‘Do Deconvolution’ button, the results shown in Fig. 6 appear. Note that, as expected from theory [1], the Maximum Likelihood regularization criterion leads to values of k and WRSS (here k= 2.7901× 10 − 3
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
and WRSS = 87.6949) lower than those obtained by means of the Discrepancy criterion, while q(k) is higher (here q(k)= 41.9112). By checking the ‘Compare solutions’ radio button in the Deconvolution Window and suitably browsing the disk, one or more previously stored solutions can be retrieved for visual comparison with the current one. For instance, in Fig. 7, we show how the program enables the user to graphically compare the solution of Fig. 6 with that of Fig. 5, previously stored in the file named ‘solution c1.mat’ (see the file list in the right of Fig. 7). Since the deconvoluted C-peptide secretion rate of Fig. 7 seems now acceptable, the ‘Confidence Intervals’ switch in the Deconvolution Window can be activated. As a consequence, after the due computations, the graphs in the Deconvolution Window are updated and the deconvoluted
75
input together with its standard intervals (9 S.D.) is displayed in the top subplot (Fig. 8). At any time, the current deconvolution problem can be quitted by pressing on the button ‘Exit’ in the Deconvolution Window and a new one can be considered from the Welcome Window. Interestingly, from the Welcome Window, one can retrieve a previously stored solution and change only a few ingredients of the problem (for example, one decides to employ a different impulse response model, e.g. a sum of two exponentials, and leave all the other information on the output and the input as they are). As apparent from the ‘Nonnegativity constraint’ radio button in the top-right of the Deconvolution Window, the program supports also the use of the non-negativity constraint if required. In the example already shown, no non-
Fig. 7. The Deconvolution Window as it appears after asking to compare the current solution (that of Fig. 6) with a previously computed solution (that of Fig. 5 in this example), archived on the disk under the name of ‘solutionc 1.mat’.
76
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
Fig. 8. The Deconvolution Window as it appears after requiring the computation of the confidence intervals for the estimate of Fig. 6.
negativity constraint violation was encountered, so the ‘Nonnegativity constraint’ option in the Deconvolution Window was never checked.
5. Conclusions The deconvolution of physiological/pharmacokinetic data is a difficult task, and plenty of methods have been developed in the literature to deal with some of its critical issues. Recently, a stochastic deconvolution approach [1] has been proposed and used in several studies [13– 16] to handle both the general and the specific difficulties. This approach, implemented by a software prototype, was successfully compared in Ref. [9] with other three software packages [3– 5] on about 50 benchmark problems covering a
wide spectrum of physiological/pharmacokinetic problems. In this paper, we presented a user-friendly computer program, WINSTODEC, that implements the stochastic deconvolution approach [1]. The major merit of WINSTODEC is its graphical user interface, which should render its use extremely easy for both physicians and biomedical engineers without specific expertise in deconvolution. Possible future developments of the WINSTODEC program include the deconvolution of time-varying systems and the estimation of nonstationary inputs (which in physiological systems can sometimes arise in response to exogenous stimuli). For the first task, since the numerical algorithms of Section 2 apply also to time-varying systems (for example, Ref. [16]), the major effort will regard the design of a new System Window
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
alternative to that of Fig. 4 (the user must be enabled to supply a time-varying kernel). As far as the estimation of non-stationary inputs is concerned, a refinement of the theory (e.g. different regularization parameters for different intervals along the input time course) can be of benefit; see Ref. [17] for some preliminary results.
Acknowledgements This work was in part supported by NIH grant RR-11095, by ICE-CNR, by MURST project ‘‘Estimation of nonassessible parameters in physiological systems’’ and by a GlaxoWellcome Research and Development grant. The C-peptide data were kindly made available by Dr K.S. Polonsky (Department of Medicine, Washington University School of Medicine, St. Louis, MO, USA).
References [1] G. De Nicolao, G. Sparacino, C. Cobelli, Nonparametric input estimation in physiological systems: problems, methods, case studies, Automatica 33 (1997) 851 – 870. [2] G. Sparacino, C. Cobelli, Deconvolution of physiological and pharmacokinetic data: comparison of algorithms on benchmark problems, in: D.A. Linkens, E. Carson (Eds.), Modelling and Control in Biomedical Systems (Including Biological Systems), Elsevier, Oxford, 1997, pp. 151 –153. [3] R. Hovorka, M.J. Chappell, K.R. Godfrey, F.N. Madden, M.K. Rouse, P.A. Soons, CODE: a deconvolution program implementing a regularisation method of deconvolution constrained to non-negative values. Description and pilot evaluation, Biopharm. Drug. Dispos. 19 (1998) 39– 53. [4] P. Veng-Pedersen, An algorithm and computer program for deconvolution in linear pharmacokinetics, J. Pharmacokin. Biopharm. 8 (1980) 463 –481.
77
[5] D. Verotta, Estimation and model selection in constrained deconvolution, Ann. Biomed. Eng. 21 (1993) 605 – 620. [6] B.R. Hunt, Biased estimation for nonparametric identification of linear systems, Math. Biosci. 10 (1971) 215 – 237. [7] D.L. Phillips, A technique for the numerical solution of certain integral equations of the first kind, J. Assoc. Comput. Mach. 9 (1962) 97 – 101. [8] A.N. Tikhonov, Solution of incorrectly formulated problems and the regularization method, Soviet Math. Dokl. 4 (1963) 1624. [9] G. Wahba, Bayesian ‘confidence intervals’ for the crossvalidated smoothing spline, J. R. Stat. Soc. Ser. B 45 (1983) 133 – 150. [10] S. Twomey, The application of numerical filtering to the solution of integral equations of the first kind encountered in indirect sensing measurements, J. Franklin Inst. 279 (1965) 95 – 109. [11] G.H. Golub, M. Heath, G. Wahba, Generalized crossvalidation as a method for choosing a good ridge parameter, Technometrics 21 (1979) 215 – 224. [12] P. Hall, D.M. Titterington, Common structure of techniques for choosing smoothing parameters in regression problems, J. R. Stat. Soc. Ser. B 49 (1987) 184 – 198. [13] G. Sparacino, C. Cobelli, A stochastic deconvolution method to reconstruct insulin secretion rate after a glucose stimulus, IEEE Trans. Biomed. Eng 43 (1996) 512 – 529. [14] G. Sparacino, P. Vicini, R. Bonadonna, P. Marraccini, M. Lehtovirta, E. Ferrannini, C. Cobelli, Removal of catheter distortion in multiple indicator dilution studies: a deconvolution-based method and case studies on glucose blood-tissue exchange, Med. Biol. Eng. Comp. 35 (1997) 337 – 342. [15] G. Sparacino, R. Bonadonna, H. Steinberg, A. Baron, C. Cobelli, Estimation of organ transport function from recirculating indicator dilution curves, Ann. Biomed. Eng. 26 (1998) 128 – 137. [16] P. Vicini, G. Sparacino, A. Caumo, C. Cobelli, Estimation of hepatic glucose release after a glucose perturbation by nonparametric stochastic deconvolution, Comp. Methods Progr. Biomed. 52 (1997) 147 – 156. [17] Pillonetto, G., Sparacino, G., Cobelli, C., Reconstruction of non-stationary biological signals via stochastic deconvolution, in: CD-ROM Proceedings of the World Congress on Medical Physics and Biomedical Engineering, 23– 28 July 2000.
WINSTODEC:
a stochastic deconvolution interactive program for physiological and pharmacokinetic systems
Giovanni Sparacino a, Gianluigi Pillonetto a, Massimo Capello a, Giuseppe De Nicolao b, Claudio Cobelli a,* a
Dipartimento di Elettronica ed Informatica, Uni6ersita` di Pado6a, Via Gradenigo 6 /A, 35131 Pado6a, Italy b Dipartimento di Informatica e Sistemistica, Uni6ersita` di Pa6ia, Via Ferrata 1, 27100 Pa6ia, Italy Received 6 June 2000; received in revised form 15 November 2000; accepted 24 November 2000
Abstract Deconvolution allows the reconstruction of non-accessible inputs (e.g. hormone secretion rate) from their causally-related measurable effects (e.g. hormone plasma concentration). Deconvolution is challenging under several aspects both general (e.g. determination of a suitable trade-off between data fit and solution smoothness in order to contrast ill-conditioning, assessment of the confidence intervals) as well as specific of physiological systems (e.g. non-uniform and infrequent data sampling). Recently, a stochastic regularization approach has been proposed and validated to handle these difficulties (De Nicolao et al., Automatica 33 (1997) 851– 870). In this paper, an interactive program, WINSTODEC, is presented to allow the clinical investigator to easily obtain the solution of a deconvolution problem by this approach. © 2002 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Software; Inverse problems; Smoothing; Biomedical systems; Endocrinology; Regularization
1. Introduction In the quantitative study of physiological and pharmacokinetic systems, deconvolution allows the reconstruction of important but non-accessi
A preliminary version of this work was presented at the Fourth IFAC Symposium on Modeling and Control of Biomedical Systems, Karlsburg, Germany, 30 March–1 April 2000. * Corresponding author. Present address: Department of Electronics and Informatics, University of Padova, Via Gradenigo 6/A, 35100 Padova, Italy. Tel.: + 39-49-8277616; Fax: + 39-49-8277699. E-mail address: [email protected] (C. Cobelli).
ble inputs (e.g. the secretion rate of a gland, the rate of appearance of a drug into the plasma after an oral administration) from samples of their causally-related measurable effects (e.g. the time course of hormone or drug in the plasma), assuming the kinetics to be linear and the impulse response to be available (e.g. from population studies or suitably designed input –output experiments). Deconvolution is well known to be challenging under several aspects both general (e.g. determination of a suitable trade-off between data fit and solution smoothness in order to contrast ill-conditioning, estimation of the confidence intervals of the deconvoluted input) as well as spe-
0169-2607/02/$ - see front matter © 2002 Elsevier Science Ireland Ltd. All rights reserved. PII: S 0 1 6 9 - 2 6 0 7 ( 0 0 ) 0 0 1 5 1 - 6
68
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
cific of physiological systems (e.g. data sampling is often non-uniform and/or infrequent, the deconvoluted input usually must be non-negative, the kinetics can be time varying). Recently, a non-parametric stochastic deconvolution method has been proposed to handle the presented difficulties [1]. In particular, thanks to the concept of virtual grid, an input profile continuous to any desired order of time derivatives can be obtained even when the available data are collected on a sparse (possibly non-uniform) time schedule. Moreover, thanks to the Bayesian framework, regularization criteria based on the Maximum Likelihood principle and confidence intervals accounting for the bias error can be derived. Non-negativity constraints can also be included. Suitable numerical algorithms permit the quick and efficient solution of deconvolution problems in a wide variety of situations, including when the system is time varying. These algorithms are based on a singular value decomposition (SVD) strategy that, by putting the deconvolution problem in diagonal form, speeds up the trial and error procedure for the tuning of the regularization parameter that controls the balance between data fit and solution smoothness. By using a prototype software, developed within the Matlab environment, the performance of this algorithm was compared in Ref. [2] with that of other three algorithms/programs from the literature [3– 5] on about 50 benchmark problems, both real and simulated, including hormone secretion problems in spontaneous as well as stimulated conditions, drug rate of absorption problems and the ‘Hunt benchmark’ [6]. This wide platform of benchmark problems covered a wide spectrum of situations: many/a few data; frequent/infrequent sampling; high/mild ill-conditioning; bad/good signal-to-noise ratio; time variance/invariance of the dynamic system. The main outcome of the comparison in Ref. [2] was that, on average, the stochastic deconvolution approach of Ref. [1] proved to be, thanks to both its firm statistical basis and numerical efficiency of the SVD strategy, more robust than the other three algorithms considered. The aim of the present paper is to introduce the reader to the interactive program WINSTODEC
(where STODEC stands for STOchastic DEConvolution), the computational core of which is basically the same as that of the prototype employed in Ref. [2]. The major merit of WINSTODEC is the new graphical user interface, which should render it extremely easy to use, despite the sophistication of theory and algorithms, for both physiologists, physicians, and pharmacokineticists as well as biomedical engineers without specific expertise in deconvolution.
2. Methods This section provides the reader with a brief description of the theory behind the WINSTODEC software. Such an overview is not meant to be complete. A full and detailed documentation of the theory is reported in Ref. [1].
2.1. The ingredients of a decon6olution problem In a linear single-input single-output system, the input u(t) and the output c(t) are linked by the convolution integral: c(t)=
&
t
g(t,~)u(~) d~
(1)
−
where g(t,~) is the kernel of the system. If the system is time invariant, then g(t,~)= g(t− ~) and g(t) is the (unit) impulse response of the system. Independently of the deconvolution approach under adoption, the two main pieces of information required to estimate the system input via deconvolution are the samples of the output and the system kernel. For instance, consider the problem of reconstructing the secretion rate of C-peptide (equal to the secretion rate of insulin, equimolarly secreted from the pancreas) in a normal man in spontaneous conditions. The ingredients are depicted in Fig. 1. The top panel displays the (noisy) samples of the C-peptide concentration in plasma collected for nearly 20 h every 10 min (therefore, the output time-series consists of n=115 samples). The bottom panel depicts the sum of two exponentials that describes the impulse response g(t) in the same subject (solid line). The impulse response
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
69
was obtained by fitting the C-peptide plasma concentration samples (diamonds) following an intravenous pulse of biosynthetic C-peptide administered in the same subject. Usually, information on the precision of the plasma concentration data is also available, e.g. measurement errors have a constant variance (known or unknown) or a constant coefficient of variation (CV) (known or unknown). In this representative example, for instance, the CV of the C-peptide plasma concentration data is 5%.
extremely sensitive to noise in the data because of ill-conditioning. Then, the idea is to establish a trade-off between fidelity to the data and smoothness of the solution, measured by a suitable penalty functional. The input estimate is represented by an N-dimensional vector uˆ, the entries of which are the samples of the input over the virtual grid chosen by the user. If y denotes the observed vector (size n)
2.2. The stochastic regularization approach
where v is the (unknown) measurement error vector affecting the n samples of c(t) in Eq. (1), and G is the n× N ‘convolution matrix’ obtained by suitably discretizing the system impulse response. Then, the estimate uˆ can be computed as
WINSTODEC implements a regularization approach [7,8]. In fact, the most direct solution of the input estimation problem (i.e. the input that, once convoluted with the impulse response, guarantees perfect adherence to the output data) is
Fig. 1. The ingredients of a representative deconvolution problem (reconstruction of C-peptide secretion rate in a normal man). Top, C-peptide plasma concentration; bottom, C-peptide impulse response (solid curve) obtained by fitting a sum of two exponentials against C-peptide samples (diamonds) following an intravenous pulse of biosynthetic C-peptide administered in the same man (for both panels, time units are in minutes).
y = Gu +v
uˆ = (GTB − 1G+ k FTF) − 1GTB − 1y
(2)
(3)
where B is an n×n matrix equal, apart from a scale factor | 2, to the covariance matrix of v, F is an N× N penalty matrix, and k\ 0 is the socalled regularization parameter. In keeping with the theory of regularization, devised in the early 1960s by Phillips [7] and Tikhonov [8], matrix F is such that Fu is the vector of the mth time differences of u, m being a positive integer parameter. In other words, F= Z m, where Z is the square n-dimension lower-triangular Toeplitz matrix, the first column of which is [1, − 1, 0, …, 0]T. In the literature, there are only a few contributions aimed at theoretically addressing the choice of m (for example, Ref. [9]). A possible reason is that the deconvoluted profile is usually not very sensitive to the value of m. As a consequence, the typical choices for m are restricted to m= 1 or m= 2, depending on the particular problem (the user finds out a posteriori which is the best value between the two). As far as the regularization parameter k is concerned, its choice is a classical topic of the deconvolution literature, where many regularization criteria have been proposed in order to avoid subjectivity in its determination. In fact, k is a crucial parameter that governs the trade-off between data fit (measured by the squared norm of the residuals vector) and solution smoothness (measured by the energy of the mth-order time
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
70
derivatives of the input). The most known regularization criteria are heuristic or deterministic, e.g. Discrepancy/Twomey [10], Generalized CrossValidation [11], Minimum Risk [12]. The stochastic approach [1] allows the derivation of Maximum Likelihood (ML) criteria, usable when the variance (or CV) of the measurement error is either known or unknown. For instance, when the measurement error variance is known, k is iteratively adjusted until uˆ TFTFuˆ =
q(k) k
(4)
where q(k) =trace[B − 1G(GTB − 1G +kFTF) − 1GT]
(5)
are the so-called equivalent degrees of freedom. The stochastic approach also allows the estimation of the confidence intervals of the deconvoluted input from the diagonal of the covariance matrix of the estimation error e= u −uˆ : cov(e) =| 2(GTB − 1G+ kFTF) − 1
(6)
This feature represents a key advantage over the other deterministic regularization approaches previously proposed in the literature; see, for example, Ref. [6] where the difficulty of handling the bias error is explicitly admitted. The treatment of the non-uniform/infrequent sampling case relies on the introduction of the virtual grid. The virtual grid is a uniform N-size grid, chosen by the user, over which the input function is reconstructed. This grid is independent of the sampling grid (the only constraint being that the virtual grid must contain the sampling grid) and can be arbitrarily fine (i.e. N n). As the virtual grid becomes finer (i.e. N ), an input profile that is asymptotically continuous up to any desired order of time derivatives is eventually obtained, even when the available data are collected on a sparse (possibly non-uniform) time schedule. The degree of continuity reflects the value of the integer parameter m selected by the user in the definition of the penalty matrix F. For instance, the typical choices of m =2 or m = 1 lead to input functions asymptotically continuous at least up to the third or first time derivatives, respectively [1].
Suitable numerical algorithms permit the quick and efficient solution of deconvolution problems in a wide variety of situations. These algorithms are based on a SVD strategy that, diagonalizing the problem, speeds up the trial and error procedure for the determination of the regularization parameter. This will enable the WINSTODEC user to quickly compute solutions for several different values of k (corresponding to different regularization criteria). Finally, the regularization method implemented in WINSTODEC can be constrained for non-negativity. In this case, the solution does not exist in closed form and an iterative algorithm, the constrained conjugate gradient method, is employed to solve the problem: min(y−Gu)TB − 1(y−Gu) +kuTFTFu u]0
(7)
3. The software This section provides the user with an overall description of the software, the use of which will be exemplified in detail in Section 4. When WINSTODEC is launched, a Welcome Window opens that asks the user to supply the ingredients of the problem (i.e. data, system impulse response, …) and the regularization settings (i.e. regularization criterion, order of the penalty functional, …). This can be achieved by retrieving a file, previously stored on the disk, or interactively by filling in a sequence of windows. The interactive modality is that of interest in this paper. For instance, in the Output Window, the user can enter the data and the measurement error statistical description; in the System Window, the user can enter the impulse response in either a parametric or non-parametric format (i.e. amplitudes and exponents of a sum of decaying exponentials or a sequence of impulse response samples); in the Input Window, the user can specify name and units of the unknown input. Once this information has been supplied, the Regularization Tools Window opens, by which the user can specify the virtual grid over which the input must be reconstructed and the order of the penalty functional to be adopted (F in Eq. (3)).
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
Once this window has been closed, deconvolution can be invoked and results are placed in the Deconvolution Window. This is a key window in the program, since it enables the user to modify the settings of the algorithm and perform deconvolution again. In particular, the user can select from a menu different regularization criteria (Maximum Likelihood, Discrepancy/Twomey and Minimum Risk, when the measurement error variance is known, and Maximum Likelihood and Generalized Cross-Validation, when the measurement error variance is unknown) or force the regularization parameter to equal a value decided by himself/herself on a subjective basis (independent of any criterion). In addition, by operating on suitable switches, the Deconvolution Window allows the user to require the introduction of the non-negativity constraint and/or the calculation of the confidence intervals of the estimated input. The outcome of the deconvolution process consists of the estimated input profile, the reconvoluted profile and the residuals, the value of the regularization parameter together with the corresponding equivalent degrees of freedom and, when requested, the confidence intervals of the estimate. All these quantities are represented within the Deconvolution Window. Graphs of the signals are displayed by default, but the user has the possibility to invoke a window where tabular values can be inspected. A typical problem with n= 115 samples and N =1150 unknown input levels like the example of this paper is solved in less than 1 min on an ordinary 450 MHz CPU/64 Mb RAM PC (an additional minute is needed to estimate the confidence intervals). From the Deconvolution Window, the user has the possibility to access to a powerful Print Menu (which enables him/her to customize panels and labels of the plots) and save the numerical results on the disk together with the settings employed. If the deconvoluted profile is not satisfactory, from the Deconvolution Window, the user can change the settings and perform deconvolution again. For instance, he/she can change the regularization criterion, modify the penalty functional, require the confidence intervals, constrain the estimate to be non-negative. When deconvolution is repeated, the computation of the solution is usually very
71
Fig. 2. The Welcome Window allows the user to enter the ‘ingredients’ of a deconvolution problem (e.g. output samples, impulse response model) either in an interactive manner or by retrieving a file previously stored on the disk.
fast (a few seconds) since, in most cases, the changes (e.g. the regularization criterion of adoption) do not affect the previously computed SVD, which represents the most computationally-demanding step of the numerical algorithm. From the Deconvolution Window, the user can also invoke menus that allow him/her to graphically compare up to ten different solutions of the problem obtained at different times (e.g. by employing different regularization settings or impulse response models or measurement error statistics description). Once closed, a problem can be later retrieved from the Welcome Window, provided that it was saved on the disk. This feature allows the user to refine in any moment the solution of a problem (e.g. by asking the program to compute the confidence intervals or by using an improved impulse response model).
4. Example In this section, the use of WINSTODEC is exemplified by solving in steps the problem of Fig. 1. Once the program is launched, the Welcome Window opens (Fig. 2), by which a new problem can be defined. Here, we consider the case in which the data are supplied interactively. In the Welcome Window of Fig. 2, the fact that any of the
72
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
labels ‘Input’, ‘System’ and ‘Output’ appears in Italic fonts reminds the user that information is still to be supplied. By double-clicking on the Output arrow, the Output Window opens (Fig. 3), and the time-series of Fig. 1 (top) can be pasted in tabular format (time, value). Labels for name (‘C-peptide concentration’) and units (‘pmol/ml’ for the y axis and ‘minutes’ for the x axis) can be supplied. The ‘Plot’ button in the Output Window (see Fig. 3) allows the user to have an x – y scatter graph of the data entered (not shown). If a threecolumn table were supplied, the third column would be interpreted as the column of the measurement error standard deviations. When the third column is not specified, as in this example, the information on the measurement error statistics (here, CV constant equal to 5%) can be entered by a suitable menu that is automatically opened by the program. This menu (not shown) proposes five choices: (a) standard deviation constant with known value to be specified; (b) stan-
Fig. 3. The Output Window allows the user to enter samples, units and uncertainty of the measured output (the representative problem of Fig. 1 is hereafter considered).
Fig. 4. The System Window allows the user to enter either a parametric (sum of exponentials model) or non-parametric description of the impulse response.
dard deviation constant but unknown; (c) standard deviation dependent on the output through a generic law (e.g. quadratic) to be specified by the user; (d) coefficient of variation constant with known value to be specified by the user; (e) coefficient of variation constant but unknown. In this example, option (d) is selected. Once the Output Window has been filled in and closed (the font of the word ‘Output’ in the Welcome Window switches from the original Italics (see Fig. 2) to Roman), one can open the Input Window by double-clicking on the input arrow of Fig. 2. By this window (not shown), name and units of the unknown input (‘C-peptide secretion rate’, ‘pmol/min’) are entered. Once the Input Window is closed, the System Window can be opened (Fig. 4). In this window, the parameters of a sum of decaying exponentials impulse response model can be supplied interactively (in this example, amplitudes 8.947×10 − 5 and 2.999× 10 − 4, exponents 0.0228 and 0.197 for the C-peptide system of Fig. 1) or loaded from a file (e.g. one can have on the disk a library of population
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
impulse response models for a number of hormones). Alternatively, one could employ a nonparametric description of the impulse response (see the non-parametric description radio button in Fig. 4) and paste, in the window that would be opened, a two-column table containing the samples of the impulse response (not shown). These impulse response samples (which can also be loaded from a file) would be linearly interpolated in order to have a continuous-time description of g(t). In both cases, the ‘Plot’ button of the System Window (see Fig. 4) allows the user to graph the impulse response just entered (not shown). Having closed the System Window, the ‘Next’ command (see Fig. 2) becomes active. One can still revise any of the ingredients of the problem by double-clicking again on its symbol in the Welcome Window, but if ‘Next’ is double-clicked,
73
the Regularization Settings Window appears (not shown), where the virtual grid period and the order of the penalty functional can be chosen. The default values proposed by the software for the virtual grid period and the penalty functional order are, respectively, the greatest common divisor of the elements of the sampling grid and 2. In this particular C-peptide secretion problem, we let the virtual grid period be equal to 1 min (this allows the estimation of a ‘quasi-continuous’ input profile since it is approximated by a piecewise constant function with 1150 levels) and m be equal to 2 (a value that was seen to work slightly better than m= 1 in this particular problem). When the Regularization Settings Window is closed, the Deconvolution Window appears, where, after having selected from the proper menu a regularization criterion (here ‘Discrepancy’), de-
Fig. 5. The Deconvolution Window as it appears after having adopted the Discrepancy regularization criterion for the problem of Fig. 1.
74
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
Fig. 6. The Deconvolution Window as it appears after having adopted the Maximum Likelihood criterion to refine the solution of Fig. 5.
convolution can be launched. The results then appear in the Deconvolution Window in a graphical fashion (Fig. 5). In particular, three subplots are presented. The top subplot displays the deconvoluted input. The middle subplot depicts the fit, i.e. the reconvolution versus the data. The bottom subplot shows the residuals standardized by the standard deviation of the measurement errors. Interestingly, by suitably clicking on any of the subplots, one can inspect a zoom view, change the labels, the curve colour, and so on. At the right of the top subplot, the program reminds the user of the settings (i.e. m value and regularization criterion) and also reports the values achieved for k, q(k) and the weighted residuals sum of squares (WRSS) (here k= 6.6996 ×10 − 3, q(k) = 35.2466, WRSS =115.007). From the Deconvolution Window, a Print menu can be invoked (not shown),
from which one can decide how to arrange the subplots before printing the results on the printer or capturing the image. For instance, one can decide to visualize only the deconvolution profile and the reconvolution versus data (not shown). From the Deconvolution Window (Fig. 5), a ‘Save’ button is also available that allows the user to store the results on the disk with an arbitrary name (e.g. ‘solutionc 1.mat’). From the Deconvolution Window, the regularization settings can be modified and deconvolution can be launched again. For example, after having changed the regularization criterion of adoption to ‘Maximum Likelihood’ and clicked on the ‘Do Deconvolution’ button, the results shown in Fig. 6 appear. Note that, as expected from theory [1], the Maximum Likelihood regularization criterion leads to values of k and WRSS (here k= 2.7901× 10 − 3
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
and WRSS = 87.6949) lower than those obtained by means of the Discrepancy criterion, while q(k) is higher (here q(k)= 41.9112). By checking the ‘Compare solutions’ radio button in the Deconvolution Window and suitably browsing the disk, one or more previously stored solutions can be retrieved for visual comparison with the current one. For instance, in Fig. 7, we show how the program enables the user to graphically compare the solution of Fig. 6 with that of Fig. 5, previously stored in the file named ‘solution c1.mat’ (see the file list in the right of Fig. 7). Since the deconvoluted C-peptide secretion rate of Fig. 7 seems now acceptable, the ‘Confidence Intervals’ switch in the Deconvolution Window can be activated. As a consequence, after the due computations, the graphs in the Deconvolution Window are updated and the deconvoluted
75
input together with its standard intervals (9 S.D.) is displayed in the top subplot (Fig. 8). At any time, the current deconvolution problem can be quitted by pressing on the button ‘Exit’ in the Deconvolution Window and a new one can be considered from the Welcome Window. Interestingly, from the Welcome Window, one can retrieve a previously stored solution and change only a few ingredients of the problem (for example, one decides to employ a different impulse response model, e.g. a sum of two exponentials, and leave all the other information on the output and the input as they are). As apparent from the ‘Nonnegativity constraint’ radio button in the top-right of the Deconvolution Window, the program supports also the use of the non-negativity constraint if required. In the example already shown, no non-
Fig. 7. The Deconvolution Window as it appears after asking to compare the current solution (that of Fig. 6) with a previously computed solution (that of Fig. 5 in this example), archived on the disk under the name of ‘solutionc 1.mat’.
76
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
Fig. 8. The Deconvolution Window as it appears after requiring the computation of the confidence intervals for the estimate of Fig. 6.
negativity constraint violation was encountered, so the ‘Nonnegativity constraint’ option in the Deconvolution Window was never checked.
5. Conclusions The deconvolution of physiological/pharmacokinetic data is a difficult task, and plenty of methods have been developed in the literature to deal with some of its critical issues. Recently, a stochastic deconvolution approach [1] has been proposed and used in several studies [13– 16] to handle both the general and the specific difficulties. This approach, implemented by a software prototype, was successfully compared in Ref. [9] with other three software packages [3– 5] on about 50 benchmark problems covering a
wide spectrum of physiological/pharmacokinetic problems. In this paper, we presented a user-friendly computer program, WINSTODEC, that implements the stochastic deconvolution approach [1]. The major merit of WINSTODEC is its graphical user interface, which should render its use extremely easy for both physicians and biomedical engineers without specific expertise in deconvolution. Possible future developments of the WINSTODEC program include the deconvolution of time-varying systems and the estimation of nonstationary inputs (which in physiological systems can sometimes arise in response to exogenous stimuli). For the first task, since the numerical algorithms of Section 2 apply also to time-varying systems (for example, Ref. [16]), the major effort will regard the design of a new System Window
G. Sparacino et al. / Computer Methods and Programs in Biomedicine 67 (2002) 67–77
alternative to that of Fig. 4 (the user must be enabled to supply a time-varying kernel). As far as the estimation of non-stationary inputs is concerned, a refinement of the theory (e.g. different regularization parameters for different intervals along the input time course) can be of benefit; see Ref. [17] for some preliminary results.
Acknowledgements This work was in part supported by NIH grant RR-11095, by ICE-CNR, by MURST project ‘‘Estimation of nonassessible parameters in physiological systems’’ and by a GlaxoWellcome Research and Development grant. The C-peptide data were kindly made available by Dr K.S. Polonsky (Department of Medicine, Washington University School of Medicine, St. Louis, MO, USA).
References [1] G. De Nicolao, G. Sparacino, C. Cobelli, Nonparametric input estimation in physiological systems: problems, methods, case studies, Automatica 33 (1997) 851 – 870. [2] G. Sparacino, C. Cobelli, Deconvolution of physiological and pharmacokinetic data: comparison of algorithms on benchmark problems, in: D.A. Linkens, E. Carson (Eds.), Modelling and Control in Biomedical Systems (Including Biological Systems), Elsevier, Oxford, 1997, pp. 151 –153. [3] R. Hovorka, M.J. Chappell, K.R. Godfrey, F.N. Madden, M.K. Rouse, P.A. Soons, CODE: a deconvolution program implementing a regularisation method of deconvolution constrained to non-negative values. Description and pilot evaluation, Biopharm. Drug. Dispos. 19 (1998) 39– 53. [4] P. Veng-Pedersen, An algorithm and computer program for deconvolution in linear pharmacokinetics, J. Pharmacokin. Biopharm. 8 (1980) 463 –481.
77
[5] D. Verotta, Estimation and model selection in constrained deconvolution, Ann. Biomed. Eng. 21 (1993) 605 – 620. [6] B.R. Hunt, Biased estimation for nonparametric identification of linear systems, Math. Biosci. 10 (1971) 215 – 237. [7] D.L. Phillips, A technique for the numerical solution of certain integral equations of the first kind, J. Assoc. Comput. Mach. 9 (1962) 97 – 101. [8] A.N. Tikhonov, Solution of incorrectly formulated problems and the regularization method, Soviet Math. Dokl. 4 (1963) 1624. [9] G. Wahba, Bayesian ‘confidence intervals’ for the crossvalidated smoothing spline, J. R. Stat. Soc. Ser. B 45 (1983) 133 – 150. [10] S. Twomey, The application of numerical filtering to the solution of integral equations of the first kind encountered in indirect sensing measurements, J. Franklin Inst. 279 (1965) 95 – 109. [11] G.H. Golub, M. Heath, G. Wahba, Generalized crossvalidation as a method for choosing a good ridge parameter, Technometrics 21 (1979) 215 – 224. [12] P. Hall, D.M. Titterington, Common structure of techniques for choosing smoothing parameters in regression problems, J. R. Stat. Soc. Ser. B 49 (1987) 184 – 198. [13] G. Sparacino, C. Cobelli, A stochastic deconvolution method to reconstruct insulin secretion rate after a glucose stimulus, IEEE Trans. Biomed. Eng 43 (1996) 512 – 529. [14] G. Sparacino, P. Vicini, R. Bonadonna, P. Marraccini, M. Lehtovirta, E. Ferrannini, C. Cobelli, Removal of catheter distortion in multiple indicator dilution studies: a deconvolution-based method and case studies on glucose blood-tissue exchange, Med. Biol. Eng. Comp. 35 (1997) 337 – 342. [15] G. Sparacino, R. Bonadonna, H. Steinberg, A. Baron, C. Cobelli, Estimation of organ transport function from recirculating indicator dilution curves, Ann. Biomed. Eng. 26 (1998) 128 – 137. [16] P. Vicini, G. Sparacino, A. Caumo, C. Cobelli, Estimation of hepatic glucose release after a glucose perturbation by nonparametric stochastic deconvolution, Comp. Methods Progr. Biomed. 52 (1997) 147 – 156. [17] Pillonetto, G., Sparacino, G., Cobelli, C., Reconstruction of non-stationary biological signals via stochastic deconvolution, in: CD-ROM Proceedings of the World Congress on Medical Physics and Biomedical Engineering, 23– 28 July 2000.