P. J. STOLL** Department of Electrical Engineering University of California Davis, California AND
J. S. MEDITCH Information Sciences Laboratory Boeing Scierdfic Research Laboratories Seattle, Washington Communicated
by Fred S. Grodins
ABSTRACT In the estimation of physiological parameters, visual fitting of experimental data has the obvious drawback that a given “best-fit curve” is not equally satisfying to every observer. In this article an iterative, weighted -onlinear least squares method of parameter estimation is formulated. It provides a s,ystematic procedure for the reduction of physiological data and obviates the need for visual fitting, which is subjective. The goodness of fit is evaluated in terms of a weighted least squares error criterion. This method is applied to estimate the parameters of a portion of the human respiratory control system. In particular, the subsystem examined is that relating tidal volume to alveolar CO8 fraction. It is modeled by a transfer function that involves five parameters : two gain constants, two time constants, and a pure time delay; and is of the same form as that determined in an earlier study involving a visual fit. The estimation is based on sinusoidal steady-state magnitude and phase data for two human subjects. The details of the numerical procedure are discussed and possible extensions are indicated. The present method improves the goodness of fit by a factor of 3 to 4. Values of the parameters changed slightly, but not insignificantly compared with those previously
* This work was supported in part by Cardiopulmonary Physiology Training Grant HEO5819-02 of the Public Health Service and in part by the Boeing Company. It was also presented at the 1lth Joint Automatic Control Conference, Atlanta, Georgia, June 25, 1970. * * Former afiliation: Division of Respiratory Diseases, Department of Medicine, University of Washington, Seattle,
opyright 0 1970 by American Elsevier Publishing Company, Inc.
P. J. STOLL AND J. S. MEDITCH
obtained using visual fitting. Implications of the results are discussed in terms of the contributions to the total ventilatory response of i,ndividual variations of CO, in brain arterial blood and cerebrospinal fluid.
Simulation of the respiratory control system has been very actively pursued [l-l I]. The resulting models h*de in common the attempt to base
a simulation of the system upon reasonable physiological and physical assumptions and upon experimentally determined subsystem parameters, whenever these are available. The inverse problem, that of estimating the system parameters from simultaneous
MODEL OF RESPIRARORY
The data presented here represent magnitude and phase of the sinusoidal component, at the driving (fundamental) frequency, of tlhe envelope described by siPLnples of breath-by-breath records of tidal volume and alveolar COz fraction. These data were obtained in earlier experiments [ 161where healthy, resting human subjects were given, for a period of lo-20 minutes, a sinusoidally varying waveform of COz at fixed amplitude and frequency. The data cover a frequency range of about 1.5 decades (0.0113 1 to 0.4582 radlsec). Details of the experimental methods for applying, recording, and estimating the magnitudes and phases of the sinusoidal envelopes of the dalta are discussed elsewhere [ 161. Magnitudes and phases of the fundamental, second, and third harmonics were all calculated, but second and third harmonic contaminations of the output waveforms were shown to be negligibly small. Thus the fundamental-frequency data were a sufficient basis for fitting a transfer function . This transfer function, denoted G(s), has the general form
where s is complex frequency, the Ki are gain constants, ri time constants, and Ti time delays. One subsystem of the respiratory regulator having a transfer function of this form is the Vr,(s)/E’Js) part of the general model shown in Fig. 1. The subsystem transfer function is denoted H,(s) in Fig. 1 and depicted in detail in Fig. 2. We then have
f-w = =
-------I71s -I- 1
(KIT2 $_ K271)S+ (Kl $_ K2) e-T3 (71s+ 1)(72s+ 1)
Obviously, H,(s) is a function of the five parameters K,, K2, rl, 72, and T.
Determination of expressions for tile magnitude and phase of H,(s) for sinusoidal frequency s = jcu is straightforward and yields the two relations
P. J. STOLL
AND J. S. MEDITCH
FIG. 1. Block diagram of CO2 regulation as determined from published data . Determination of hypothetical transfer relationship W from F4(s) to d(s) was beyond the scope of the reference and is conjectural at this time. Of particular interest in this report is the V,‘,B(~)/F,(~)transfer function HI(s). A full list of symbols appearing in this diagram is given in Ref .
(K,2+ppw_ fan-l 1
(71+ 72bJ- WT, 1
respectively. The experimental data from which the parameters will be estimated are given in Table I and were obtained in the following way [ 161. The experiment was repeated several times at each test frequency for each of the two subjects studied. Each magnitude or phase value at a given frequency was regarded as a datum. The mean and variance of the collection of magnitude (phase) values at each frequency were then computed for each subject. In the previous study , a smooth curve was visually fitted to the mean values for each subject as indicated in Fig. 3, and the parameter Mar hema tical Biosciences
Fluid and Tissue
FIG. 2. Detailed structure of V”,(s)/F’(s) transfer function H,(s). Short time constant q reffects mainly changes in arterial CO%concentration, while the longer time constant 7a probably reflects some composite effects of much slower changes in CO, concentration in brain tissue and cerebrospinal fluid. TABLE I RECORDED MAGNITUDE
AND PHASE DATA
AND NUMBER OF REPETITIONS (n)
FOR SUBJECTS 1 AND
0.01131 0.01707 0.02775 0.04370 0.07016 0.1121 0.1791 0.2853 0.4582
0.5252 0.5044 0.4350 0.3710 0.2585 0.1646 0.1155 0.08658 0.1380
A.6248 -0.7978 -0.9192 -1.165 -1.511 -1.969 -2.747 -3.269 -4.060
0.006247 0.01329 0.01968 0.007669 0.003507 0.006250 0.003787 0.003898 0.005201
0.006084 0.006849 0.004324 0.01716 0.01266 0.02396 0.7585 2.028 1.145
6 7 8 7 8 8 7 8 7
0‘01131 0.01707 0.02775 0.04370 0.07016 0.1121 0.1791 0.2853 0.4582
0.6924 0.7424 0.6903 0.6372 0.5099 0.2891 0.2121 0.1105 0.1511
-0.4189 -0.6135 -0.6250 -0.8674 -1.199 - 1.740 -2.176 -3.063 -3.227
0.03889 0.01633 0.06520 0.03437 0.003864 0.004981 0.003166 0.005580 0.006964
0.03119 0.01988 0.04610 0.01182 0.01748 0.04700 0.4477 4.219 0.5616
6 6 6 6 6 6 6 5 3
Subject 1 1 2 3 4 5 6 1 8 9 Subject 2 1 2 3 4 S 6 7 8 9
Mathematical Biosciences 8 (1970), 307-321
P. J. STOLL AND J. S. MEDITCH
0 6 CL
0, ‘; d
003 002 001
I I 03
1 / -50
-.-. . / 01
j, 32 Podtans
FIG. 3. Frequency response diagrams f’or I/TA(s)&(s). A, magni:dde and phase for ect 1; B, magnitude and phase for subject 2. Solid curves are fitted visually to mean es of data and are plots of E.qs. (3) and (4) with values of parameters given on first line (i = 0) of Tables II and IV, respectively. Dashed curves represent the same transfer function [Eqs. (3) and (4)] plotted with estimated parameters from Table II, line 4, for subject f (Fig. 3, left) and Table IV, line 7, for subject 2 (Fig. 3, right).
estimates were determined as follows. First, a single time delay T was obtained to correspond with the high-frequency portion of the phase versus frequency data, together with a single time constant and gain to fit e magnitude versus frequency curve. This procedure utilizes the fact at the transfer function was minimum phase apart from the time delay. second time constant and gain were added to improve the fit, but addiof more time constants and gains resulted in little fuxher improvet. The two gain constants KI and KS were then determined from the low-frequency gain is KI + Kz. re utilizes only the mean values of
the measured magnitudes arid phases. No relative weighting was attached to these quantities in terms of their variances, which are indicative of the confidence associated with each such mean value. In the present study, each mean value will be treated as a measurement and its corresponding variance as descriptive of the error associated with that measurement. The variances will then be used as weights in the least squares estimation performance criterion. Also, we will use the parameter estimates obtained previously  as the “starting values” for the iterative estimation procedure described in the next section. These estimates are given in the first lines of Tables II and IV for subjects 1 and 2, respectively. PARAMETER
The problem of determining the set of parameters (&, K,, 71, a2, T) for each of the aforementioned two subjects can readily be formulated as one of nonlinear weighted least squares estimation [ 171. To facilitate the development, we consider a general problem of estimating a set of parameters x where x is an n-vector that is nonlinearly related to a set of measurements
z(w3 = h(x, mk)+ ebk)
wherek = 1,. . . , p. Here, ~r)~is the discrete independent variable, usually time, but frequency in our case; h is an m-dimensional vector-valued function of the indicated variables; z is the observed measurement, an m-vector; e is the measurement error, also an m-vector; and p denotes the number of points at which measurements are given. The nonlinear least squares estimation problem consists in determining x such that the quadratic form VR =
- W, qJ1
is minimized. In Eq. (6), the prime denotes the transpose, [ 1-l the matrix inverse, and R(o,), k = 1, . . o , p, is a set of m x ZVE weighting matrices. These matrices are selected to refrect one’s confidence in each measurement. Clearly, T/n is a scalar function of the residuals z(cuk)
and is a quadratic measure of the extent to which an estimate of x fits the observed data.
P. J. STOLL AND J. S. MEDITCH
For our particular case, n = 5, m ==2, and p = 9 with .
*1 x2 x-c
Hence, from Eqs. (3) and (4), we have (Xl
M*, qi.) = IWco,)l =
+ I +
63 _. + x,)Wk 9L - x3xp;
Xl + x2
where k = 1, . . . ,9.
For the weighting matrices, we have chosen
9 where ~‘&(~~,J and c$,(wJ are given in Table I for the two cases of interest. Since kl(C(Jk) appears in vR in Eq. (6), each term in VR is weighted according to the reciprocals of the two variances, 0’41(%)and $,(%)* Finally, since zl(cok)and z&+) are the measured values of lH’(jc+)l and /Hl(j&, respectively, we have zl(cr)k) = hl (x,
and z2(Wk)= h&9 %) +
where h,(x, 0~~)and h,(x, wk)are defined above, and &&.) and &ok) are the corresponding measurement errors. Returning to the general problem, we set the gradient of VR with respect to x equal to zero to obtain the result *, k=l
6Jk)R’%Jk)[Z(6Jk) - h(x,
a,)] = 0.
n this relation, hJx, 6~~)is the m x n Jacobian matrix of h(x, (ok)with
RESPIRATORY SYSTEM PARAMETERS
respect to X; its elements are ah,& w&/ax, where A = 1, . . . , m denotes the row of the element and a = 1, . . . , n the column. Equation (9) is a system of n nonlinear algebraic equations in as many unknowns and is a necessary con,dition for x to minimize VB. If we dencte this system of equations by gO4 = 0, (lo) then Eq. (9) is also a sufficient condition for vB to be a minimum if the n x n matrix g.(R) is positive definite, where R is a solution of Eq. (10). Straightforward calculation gives
where B(x) is an n x n matrix whose elements zre given by the relation
whereaJ= l,..., n, and a and /I denote the row and column, respectively, of the element. Also [R-‘(&)]~, is the ilq element of R-‘(ok) and [z(wk) - h(x, cl)k)]qiSthe qth COmB,Onent Ofthe wt-Vectorz(mk) - h(x, ok). We note that g,(x) is a symmetrical matrix. Since the system of equations in Eo, (lo), or equivalently Eq. (9), is nonlinear, the system cannot be solved in closed form except in certain special cases. It is, therefore, logical to utilize Newton’s method to formulate an iterative procedure for determining the required roots 1181. Letting the iteration index i = 0, 1, . . . , Iv, where N is some integer, we assume that the ith iterate xi is given. Then, expanding the left-hand side of Eq. (10) in a Taylor series about xi and retaining only first-order terms, we get the linear system of equations
g(x”) + for the (i + 1)th iterate.
xi] = 0
If the process defined by Eq. (13) converges at all, it will do so quadratically [ 181. However, the conditions that must be satisfied to ensure convergence are virtually impossible to check in practical problems of even modest dimension. We note that the iteration process involves both first and second partial derivatives of the components of h with respect to the elements of x. Thus, the method is second order in measurement function nonlinearity. This is in contrast to the usual procedures, which involve a priori linearization , and are, therefore, only first order. The iteration process is Mathematical Biosciences
P. J. STQLL AND J. S. MEDITCH
started with an “initial guess” x0 and may or may not converge, depending on the particular value chosen. A digital computer program was developed to implement Eq. (13)
for estimating the parameters (K, K,, To, TV, T) using the data in Table I. Numerical differentiation using &e-point difference approximations was employed to compute the partial derivatives a&(x, w,)/ax, and a2hA(x,cok)~i3x,ikfl,a, ,8 = 1, .
. . ,
5, a = 1,2,
needed in the evaluation of g(x”) and gz(xi). This introduced noise into the computational process as shown in Tables II-IV. Further, a certain amount of trial and error was necessary in the choice of difference increments in order to obtain reasonable behavior in the computations. computational results for the measurement data of subject 1, Table I, are given later in Tables II and III. In Table II, previously determined parameter values 1161were used for the initial iterat , whereas in Table III, the parameter values for the initial iterate were selected arbitrarily, but in a neighborhood of the values to which the parameter estimates converged in Table II. The value of YR, as defined in Eq. (6) is included in these tables to show the magnitude of the estimation performance criterion as a function of the iteration index. All results are presented to four significant figures, since all input data were given to the same accuracy. The calculations were:carried out on an IBM 360/44. Double-precision arithmetic was used to irvoid the possibility of roundoft’ error accumulation. However, the e%ct of error introduced by using single-precision arithmetic was not e::?Gned. 1Jsing the behavior of VR as a function of the iteration index d as the
criterion of convergence, the process is seen to converge in four iterations in Table 11 and three iterations in Table III to within the indicated accuracy. A pl’ot of the transfer function using the corresponding parameter values is included in Fit:;. 3A. Allowing the computations to continue to ten iterations clearly poiints out the noise introduced by numerical differentiation, the effect being most pronounced in the estimation of Kz and 72. This relative insensitivity of vR with respect to K2 and TVwas characteristic of all of our calculations. This is probably due either to inaccuracy because of the finite difference approximation of the partial derivatives or to the dominance of the transfer function by K,, Q, ard T. Ineither event, the minimum value of Vn appears “sharp” in &, TV, and T, but “flat” in Kk and r2. Table :W gives the computed results for subject 2. Here, vn varies between 3.195 and 3.197 during the last five iterations, The corresponding
uctuations in K2 and r2 are more marked than those in the other paramst as for subject igure 3 lot of the transfer .Wathema!ical Biosciences
FUNCTIONS OF IT’BRATION INDEX i FOR SUBJECT 1’
0 1 2
0.4419 0.447 1 0.4473
3 4 5
3.024 3.024 3.024
a Initial estimate (i = 0) from previous work .
function with the estimated parameters corresponding to minimum VB. In Table IV, the value of vR actually increases for the first two iterations before recovering. This indicates that Newton’s method “overcorrects” initially, a phenomenon most likely avoided by making
FUNCTIONS OF ITERATION INDEX
& (liters/ %CO,)
& (liters/ %CO,)
0.4000 0.439 1
6.21 I 6.211
0.447 1 0.4461
o Initial estimate arbitrary;
see text. Mathematical Biosciences
J. P. STOLJ, AND
J. S. MEDITCH
& (liters/ WACO&
26.60 3.949 9.925
u Initial estimate from previous work .
4 for example. This forces Newton’s method
The noise introduced into our computational results as a consequence of numerical differentiation to evaluate g(x”) and g,(x”) has been noted. In retrospect, we would have perhaps encountered less difficulty if the required partial derivatives had been determined analytically. In our case where n = 5 and HZ= 2, this involves 40 partial derivatives of which 11 are identically zero because of the functional dependence of !zl and kz on x. Although tedious, the task is feasible. However, for a more complex model, for example where n = 12 (i.e., 12 parameters) and 172= 2 again, 180 partial derivatives are required. Unless the functions hl and ha are particularly simple and a substantial percentage of the required partial derivatives can be shown to be identically zero, we are forced to use numerical differentiation to make the problem tractable. No attempt was made here to optimize the numerical differentia”ion process. For example, we used the same increment sizes to compute both the first and second partial derivatives. Further attention to this aspect of computation would lead to more desirable results, not only in the present case, ut also in future applications. (1970), 307-321
RESPIRATORY SYSTEM PARAMETERS
In a similar vein, certain variants of Newton’s method, termed “extended Newton’s methods” [2O], offer some promise of improving the computational procedure. The simplest such modification or extension has been noted in Eq. (14), for which a number of criteria for selecting p are possible. Since Eq. (9), or equivalently Eq. (lo), is nonlinear, uniqueness of the solution obtained by Newton’s method, or any other method, cannot be guaranteed. The question still remains whether or not the solution is a relative minimum of VR. In such cases, we must rely upon physical intuition and insight into the particular problem at hand. The presence of the pure time delay means that V,(S)/F~(S) is a nonminimum-phase transfer function. This, in turn, implies that phase as well as magnitude data are necessary to estimate the parameters. On the other hand, a minimum-phase transfer function would be uniquely determined by its gain characteristic alone. Earlier in the text and in Fig. 2 the V,,(s)/FJs) transfer function was seen to consist of a “fast” component (K,, Q) corresponding to arterial COZ in the medulla oblongata, and a “slow” component (kl,, Q) for cerebrospinal fluid (CSF) CO2 variations. This reasoning is based on work by Pappenheimer and others , who argued that ventilation is a single function of the hydrogen ion concentraticn at a location within the brain tissue where the effect of changes in hydrogen ion concentration is intermediate between variations of this ion in brain arterial blood and CSF. During COZ breathing, changes in hydrogen ion concentration, to which the brain chemoreceptors respond, rapidly follow COZ changes. In the model presented here, the single time delay (T) represents the transport time for blood to move from the lungs to the medulla. The time constants Q and 72 correspond to the fast and slow rates of buildup of COz and hydrogen ion concentration at receptor sites exposed primarily to arterial blood and the less accessible CSF, respectively. Other functional forms of the V&)/J”‘(s) transfer function could be tried and the corresponding magnitudes of Vn compared with those obtained here. The best-fit functional form and set of parameters would yield the minimum yIz and yet be reasonable physiologically. As yet, these forms would be selected on a trial-and-error basis. Remaining to be developed is a means for simultaneously optimizing both the mathematical form and values of the parameters of the transfer function that best fits the data. One possible function for further trials was obtained from step-response data by Lambertsen et al. [ 131and contains a term [K3exp( - T2s)]/(~3~-I- 1) in addition to those in Eq. (2~. The time constant 73 and time delay IT2are both very short, and corres onse of the arterial chemoreceptors Mathematical Bio.viences 8 (1970), 307-321
J. P. STOLL
AND J. S. MEDITCH
to a sudden change in arterial COz. This term contributed much less to the total ventilatory response than the remainder of the transfer function, amounting to only 12% of the total; it was not visually discernible from the sinusoidal data and therefore not included in either the visual or the least squares estimation processes given here. During the parameter estimaticn, the values of &, Q, and Q changed but slightly from the initial estlnrates for the first subject (Table II); & and T changed considerably rnkjre. For subject 2 (Table IV), Kg and Q changed slightly, while changes in the other parameters were relatively larger. The near-zero frequency gain was increased from about 1.11 to ZA liters/ %COZfor subject 1 and &creased from 2.07 to 1.68 for subject 2 during the iterative process. The estimated curves still do not differ greatly from the visually fitted curves (Fig. 3) within I’llerange of experimental data; the discrepancy between dc gains results from the different behavior of the fitted functions at frequencies below the range of data (below 0.01 rad/sec). Until further data in the 0 to 0.01 rad/sec range are obtained, significance cannot be assigned to the disc :pancy in dc gains. We thus see that the fit most pleasing visually to an observer can be close to, but not necessarily identical with, the one giving minimum error. When the transfer functions were visually fitted to the data, only t overall dc gain and time delay differed between subjects 1 and 2; K&K& rl, and 72 were the same for botk subjects. The conclusion was thus reached earlier that the same form of transfer relationship described tidal volume response in both subjects, with differences only in the circulatory delay (T) and overall sensitivity. With parameter estimation, however, all parameters differ between the two subjects, although the functional forms are the same. In subject 1, = 0.463 and in subject 2, 0.738. This compares with K,/& of WK, 0.795 for both subjects in the initial estimates, The parameter estimation method resulted, therefore, in more weight being given to the longer time constant portion (acting through CSF) of the system relationship. However, what differences exist do not change the gain-phase plots greatly. In view of these small dif%rences between the time constants for the two subjects, we can still safely say that the same mechanism is operating in both cases; individual variations in physiological constants between different subjects can account for the differences between the fitted parameters.
1 J. G. Det’ares, X-I.E. Derksen, and J. W. Duyff, Cerebral bloodflow in the regulation of respiration, Acta Physiol. Pharmacol. Neerl. 9( 1960), 327-360.
RESPIRATORY Sk STEM PARAMETERS
2 F. S. Grodins, June Buell, and A. J. Bart, Mathematical analysis and digital simulation of the respiratory control system, J. App/. Physiof. 22(1967), 260-276. 3 F. S. Grodins, J. S. Gray, K. R. Schroeder, A. L. Norins, and R. W. Jones, Respiratory responses to COB inhalation. A theoretical study of a non-linear bioor,J, Appl. Physiol. 7(1954), 283-308. 4 F. S. Grodins and G. James, Mathematical models of respiratory regulation, Ann. N. Y. Acud. Sci. 109(1963),852-868. 5 A. C. Guyton, J. W. Crowell, and J. W. Moore, Basic oscillating mechanism of Cheyne-Stokes breathing, Amer. J. Physiol. 187(1956), 395-398. 6 J. D. Horgan and R. L. Lange, Digital computer si ulation of respiratory response to cerebrospinal fluid PC0 in the cat, Biophys. J. 5 7 R. L. Lange, J. D. Horg!an, J. T. Botticelli, T.
Kuida, Pulmonary to arterial circulatory transfer function: Importance in respiratory control, J. Appl. Physiol. 21(1966), 1281-1291. 8 H. T. Milhorn, R. Benton, R. Ross, and A. C. Guyton, A mathematical model of the human respiratory control system, Biophys.J. 5(1965), 27-46. 9 H. T. Milhom and A. C. Guyton, An analog computer analysis of Cheyne-Stokes breathing, J. Appl. Physiol. 20(1965), 328-333. 10 C. M. E. Matthews, G. Laszlo, E. J. M. Campbell, P. M. Kibby, and S. Freedman, Exchange of “CO8 in arterial blood with body CO2 pools, Respir. Physiol. 6(1968), 29-44.
11 C. M. E. Matthews, G. Laszlo, E. J. M. Campbell, and D. J. C. Read, A model for the distribution and transport of CO%in the body and the ventilatory response to CO:!, Respir. Physiol. 6(1968), 45-87. 12 P. Padget, The respiratory response to carbon dioxide, Amer. J. Physiol. 83(1928), 384-394.
13 C. J. Lambertsen, R. Gelfand, and R. A. Kemp, Dynamic response characteristics of several CO,-reactive components of the respiratory control system, in Cerebrospitlalfluid and the regulation of veratilation(C. MC. Brooks, F. F. Kao and B. B. Lloyd, eds.), pp. 211-240. F. A. Davis, Philadelphia, 1965. 14 E. Florentin, Etude theorique et analogique de la composition du gaz alveolaire en regime periodique. Application au probleme de la regulation chimique de la respiration, J. Physiof. (Paris) 56(1964) (Supple. IX), l-60. IS J. P. Thompson, Alralysisof respiratory control using sinusoidally varying inspired carbon dioxide. Master’s thesis, Univ. of Washington, Seattle, 1962. 16 P. J. Stoll, Respiratory system analysis based on sinusoidal variations of COz in inspired air, J. Appl. Physiof. 27(1969), 389-399. 17 K. F. Gauss, Theory of the motion of the heavenly bodies moving about the sun in conic section? (reprint ; C. H. Davis, trans.) Dover, New York, 1963. 1q A. A. Goldstein, Constructive real analysis. Harper & Row, New York, 1967. 19 R. Deutsch, Estimation theory. Prentice-Hall, Englewood Cliffs, New Jersey, 1965. 20 L. Collatz, Functional analysis and numerical mathematics, (H. Oser, trans.). Academic Press, New York, 1966. 21 J. R. Pappenheimer, V. Fencl, S. R. Heisey, and D. Held, Role of cerebral fluids in control of respiration as studied in unanesthetized goats, Amer. J. Physiol. 208(1965), 436-450. Mathematical Biosciences