1 Introduction THE DYNAMIC ventilatory response to hypercapnia and exercise has been investigated by frequency responses (DAUBENSPECK,1973; CASABURIet al., 1977; SWANSONand BELLVILLE, 1974; WIGERTZ, 1970) and impulse responses (BENNETT et al., 1981; SOHAB and YAMASHIRO,1980) as well as transient responses to various inputs. However, the dynamic response to their combination has not been analysed. Moreover, previous investigators seemed to pay relatively little attention to inherent system noise. (An exception is the issue of how it can be eliminated from the estimated ventilatory response.) As ventilation shows considerably irregular fluctuation, it is important to know the characteristics of noise as well as the system (plant and controller) for the analysis of dynamic ventilatory response. It is especially crucial in an attempt to predict and control the subject's ventilation. To accomplish this, it is necessary to decompose ventilatory fluctuation into the inherent noise and influences from other variables. Multivariate autoregressive modelling and relative power contribution analysis are useful methods of discriminating these components. The relative power contribution analysis decomposes the power spectrum of ventilation into the sum of the contributions from each variable at each freCorrespondence should be addressed to Dr Yoshitaka Oku, Department of Biomedical Engineering, Case Western Reserve University, Cleveland, OH 44106, USA. First received 17th August 1990 and in final form 8th March 1991

9 IFMBE: 1992

Medical & Biological Engineering & Computing

quency. The mathematical procedure is described in the Appendix. Decomposition is feasible even if the system includes multiple variables which interact with each other. This is an advantage over the cross-spectral method which gives a biased result for the analysis of the feedback system (AKAIKE,1967). In the present paper the dynamic influences of end-tidal CO2 and exercise on ventilation are examined when their levels are perturbed separately and simultaneously. Work rate and end-tidal CO2 gas concentration were changed quasirandomly, and the contributions of end-tidal CO2 concentration (FETcO2) and work rate (WR) on ventilation (l/e) were examined by the autoregressive model with the exogenous variable (ARX model). The power spectrum of ventilation was decomposed into the contributions from each variable and the residual power, which cannot be accounted for by the other variables. 2 Methods The subjects were five healthy young adults ranging in age from 24 to 29 years. The subjects were familiar with the breathing apparatus, but were not informed about the experimental protocol beforehand. While in the supine position on a bicycle ergometer (Tatebe Seishudo EM-401), they placed their feet in the stirrups on the pedals and were instructed to pedal constantly at 60 rotation min-1. Inspired and expired flow was measured by a hot-wire flowmeter. Gas fractions in the mouth were measured with a mass spectrometer (Perkin-Elmer MGA 1100). A low

January 1992

51

dead-space oronasal mask, through which room air was breathed, was fitted on the face. The total dead space, including the breathing valve and the mask, was 30-35 ml. Minute ventilation (1?E), 0 2 uptake (17o2), CO2 output (17co2), heart rate and FE~co~ were computed at intervals of 10s by a respiromonitor (Minato RM300) using breath-bybreath estimations. The study consisted of three trials: (a) random work rate forcing; (b) random CO2 inhalation; (c) random work rate forcing with random CO2 inhalation. Each subject performed only one trial per day, and the sequence of these trials was randomised among the subjects. (a) Random work rate forcing: The work load was switched over between 20 and 80W for all subjects. It was determined that this value did not exceed the anaerobic thresholds, which were estimated by a 2 min incremental work test in a manner similar to that described by WASSERMAN et al. (1973). The work load was maintained below this value because the ventilatory response for the work load is nonlinear beyond Table 1 Physical characteristics and anaerobic thresholds for each subject

Subject E.K. O.N. K.H. N.M. H.N.

Physical characteristics Age (years) Weight (kg) Height (cm) 25 65 170 27 54 171 25 73 170 29 56 171 27 65 174

AT (W) > 100 100 140 120 120

the anaerobic threshold. Anaerobic thresholds and physical characteristics for each subject are presented in Table 1. Subjects rested in the supine position with the oronasal mask present for at least 15 min, followed by a 5 min warm-up period at 20W. The work rate was then varied in a random sequence for 30min. The duration of each work load was determined quasirandomly to lie between 10 and 90 s at 10 s intervals. (b) Random COz inhalation: Inhalation gas mixtures were prepared by a random gas generator (Respy Laboratory GASMIXER-002). Pure gas sources of N 2, 0 2 and CO2 were connected to a gas mixing chamber through electromagnetic valves which were controlled by an 8-bit digital computer (Fujitsu FM-77). The system was capable of mixing gases to any desired final composition with a first-order delay of about 3 s. The to subject from subject I

I

mixed gas

I random gas generator

t

gases were humidified at room temperature as they flowed to the chamber. The concentration of inspired CO 2 was varied randomly for 30rain between 0 and 7 percent at random intervals from 10 to 90 s. The amplitudes of work rate forcing and inspired CO2 fluctuation were adjusted so that each produced approximately the same amount of ventilatory fluctuation. This CO 2 random signal was generated so that there was no cross-correlation with the pseudorandom work rate sequence. (c) Random work rate forcing with random CO2 inhalation: The same pseudorandom sequences as described above were imposed on each subject simultaneously. Fig. 1 shows the system identification process we used. The data obtained from the respiromonitor were analysed by a digital computer system (DEC VAX 11/750). The data obtained during the first 5 min of each trial were discarded and the data obtained during the subsequent 25 min. were analysed. The TIMSAC (time series analysis and control) computer program package written by AKAIKE and NAKAGAWA (1972) was used. The frequency domain characteristics of work rate fluctuation are presented in Fig. 2. This figure shows the power spectral density for work rate fluctuation. The point at which the power dropped below the 3 dB level, i.e. the half-power point, was considered to represent the upper limit of useful frequencies. This point lies approximately at the frequency of 1 cycle min-1. The power band is slightly narrower than that of the pseudorandom sequence described by BENNETTet al. (1981). Analysis of the correlation coefficients of the residuals (innovation) is useful for the selection of system variables to be used eventually. Initially, we took ventilation (1?e), gas exchange variables (Vo2, Vco2 and FETCO2) , heart rate and work rate. Heart rate was chosen to examine the influence of cardiodynamic changes on ventilation. The crosscorrelation coefficients of the residuals normalised by the autocorrelation of the residuals were as shown in Table 2. In all subjects, the residuals of 1?E, 17o2and 17co2showed a strong correlation. Several laboratories have taken the good correlation between 1?E and 17co2 as suggestive evidence of some underlying CO2 linked control mechanism that can operate independently of the mean level of PaCO 2 . However, causality cannot be discriminated between them, and it is also possible that some of the apparent correlation is due to a common signal that drives 1.0

I

random work load

I

hot wire I flowmeter

maSS spectrometer

bic r

ergometer o o.

data acquisition

0"5

1

L quasirandom signals

9 . . [ VE J VO2 JVCO2 J FETC02. heart rate ARX model fitting

l

0

i

I

I

i

2

3

frequency, cycle rain-I

relative power contribution analysis Fig. 1 Block diagram of system identification process 52

I

Fig. 2

Normalised power spectral density for the input sional (random binary work rate forcino). Upper limit of useful frequency content is that frequency at which the power falls below the 3 dB level

Medical & Biological Engineering & Computing

January 1992

I/E and l?co~. In such a case, considerable errors are introduced to the estimation of the relative power contribution, and so their simultaneous inclusion in the model must be avoided. When I?E, FErco~, heart rate and work rate were selected as variables, the cross-correlation of the residuals was small enough for practical use as shown in Table 3. Table 2 Normalised correlation coefficients between innovations (for subject N.M.)

176 .1?o2 Vco~ PExco~ HR WR

17o2 0'87 1.0

1.0

1?coz 0"96 0-94 t.0

P~TCO2 -0.27 - 0-06 -0-20 l'O

1?E= minute ventilation; I?o2, = 0 production; PETCO2= end-tidal HR = heart rate; WR = work rate

2

HR 0-25 0.42 0-31 O"13 1"0

WR 0.002 - 0-02 -0.0006 O"10 0"23 1"0

consumption; 17co2= C O 2 CO2 partial pressure;

Table3 Normalised correlation coefficients between innovations (for subject N.M.) 1?6 1-0

I26 PETCOz HR WR

HR 0'25 0"31 1"0

PExco: --0"39 l "0

WR -0-01 -- 0"04 0'23 1"0

I?6 = minute ventilation; PETCO2= end-tidal CO 2 partial pressure; HR = heart rate; WR = work rate 3 Results The raw data for subject N. M. (quasi-random work rate forcing only) are shown in Fig. 3. Frequency responses of ventilation to work load during exercise alone and simultaneous CO 2 and exercise loading are shown in Fig. 4. As the frequency of the fluctuation of work rate increased, the phase angle for ventilation to the work rate and the amplitude of the ventilatory response diminished (low-pass filtering property). These findings are similar to those of WIGERTZ (1977) and CASABURI et al. (1977). The coherency of ventilatory response to work rate greatly decreased with simultaneous CO2 and exercise loading. The relative power contribution of each variable to ventilation was given by decomposing the power spectrum of ventilation into the sum of the contributions from each variable at each frequency. Note that it is expressed as a percentage of the entire power of ventilatory fluctuation rather than an absolute value. Weighted ensemble means 80

work rate,W 2 9 4 VE , I min -I

0 0

~ ~

I0

1500r {/o2.mi min -I 0

15oor 0 ~

VC02 ml rnin-I FETC02' per cent

~

o

~

~

~

i

beats/rnin ~1

Fig. 3

time, min Changes in ventilation and gas exchange variables in a sample run during random work rate forcing

Medical & Biological Engineering & Computing

of the relative power contribution of ventilation during work rate forcing, CO2 inhalation and both are shown in Fig. 5. Each area represents the contribution of each variable to ventilation at each frequency. All subjects exhibited similar contribution characteristics. The bottom zone of each frame (12~) shows the part of the spectrum which cannot be accounted for by other variables, i.e. the contribution of the residual power. At low frequencies, this is approximately 30 per cent during either CO 2 or exercise perturbation, and rises sharply with the frequency. This indicates that the contribution of noise or nonlinearities is quite high. It does not mean that the absolute power of noise components grows as the frequency increases because the power spectrum of ventilation decays as the frequency rises. The residual power was smallest in the simultaneous CO 2 and exercise perturbation. In the case of work rate forcing alone, the contribution from work rate to ventilation was significant in the low frequency band, but decreased gradually as the frequency increased. In the case of CO 2 inhalation alone, the contribution from FETCO2 to ventilation was significant in the low frequency band and gradually decreased as the frequency increased. The relative contribution of F~xco ~ during CO 2 inhalation alone showed only a slightly larger contribution than that of the work rate during random binary work rate forcing alone. The ratio of these contributions is about 7:5 in the low frequency band. Thus, if the system is linear, the principle of superposition is applicable and the relative contributions of hypercapnia and exercise during simultaneous loading would be the same as those predicted by this ratio. However, in the case of simultaneous loading of CO 2 and exercise, only the contribution from FExco2 was significant. The relative contribution of the work rate was smaller than expected, which shows the nonlinear system dynamics for the combination of exercise and CO 2 inhalation forcing. 4 Discussion Our results have shown that the contribution of noise or nonlinearities is quite high even at low frequencies for either CO2 or exercise perturbation and that it rises sharply as the frequency increases. The residual power was smallest in the simultaneous CO2 and exercise perturbation, which would be a reflection of an increased signal-tonoise ratio. We have also shown that the relative contribution of the work rate during simultaneous loading of exercise and CO2 is smaller than the expected value when we assume that the system is linear. On the other hand, the contribution of CO2 is greater than expected. This observation means that the combination of exercise and CO z forcing reveals nolinear system dynamics, and ventilation tends to follow the fluctuation of CO2 level when the work rate and the CO2 level vary simultaneously. It is consistent with the decrease in coherency observed when exercise is combined with CO2 inhalation (Fig. 4) which indicates an increase in noise and/or nonlinear behaviour. A pseudorandom test signal is useful for the dynamic system identification and several investigators applied this to the respiratory system (BENNETT et al., 1981; SOHAB and YAMASHIRO, 1980). An application to the cardiovascular system has also been reported (SUYAMAe t al., 1988). Random forcing of inhaled CO 2 concentration gives enough power over a wide frequency range. We have not imposed a quasirandom signal on the end-tidal Pco2, but rather on the inspired Pco2. Thus, the end-tidal signal is more frequency-dependent. However, the present method does not require whiteness of the input variables. The only prerequisite is that the residuals are uncorrelated. This is

January 1992

53

0

160

~

ao

u>, 0-6

"0

-E 2O ._~ 10

5 0

0.8

-60 ~ [ ~

T ._~ E 40

~

cycle min-! 0'5 1-0

r

-120

c-

-180

~, O'a 8

t3.

-24q

0-5 cycle rnin-1

0.2 0

12

cycle min-I 0.5

320

0.'5~0 cycle rnin-I

1.0 0.8

.o -60

0.6 -120

_~ 4o

o 0.4 o

.~ 20

~, -~a0 c-

1-o

50

Fig. 4

0.2

O.

0'5

-24(

I~0

1.~0.50 cycle min-I cycle min-I Weighted ensemble means of gains (amplitude ratio of input and output), phase shifts and coherencies in frequency responses of ventilation to work load plotted in a Bode diagram. Responses were normalised on a per W basis: (a) without CO 2 inhalation; (b ) with CO 2 inhalation.

100

'::::::::i:iii!ii:.~i~iiiii::::::iii::i;:iii::iiiiiii iii~iiiiiiiii::ii~ii::ii~iiii~i~:i:~:i::::: ~:~:~:~:~:~:~:~:~:::i...............

.;,:.:o~, ,_~, 9 9 ~

..==o,1.o=1=%%* ==%%%,===,.ol%eo-e~o

t-

t). O

50 u

P 0"2 0"4 0"6 O'B 1"0 0.4 0"6 0"8 1.0 0 frequency, cycle rain-! frequency, cycle rain-I c b Weighted ensemble means of the relative power contribution for ventilation. (a) during work rate forcing without C O 2 inhalation; ( b) during C02 inhalation only; ( c ) during work rate forcing with CO 2 inhalation

0.2

Fig. 5

0-4 0-6 0.8 frequency, cycle rain-1

1-0

0

0.2

an advantage over the direct impulse response estimation in the frequency domain. As there was no cross-correlation between quasirandom signals of work rate and CO2, cross-correlation between residuals of two variables (FETcO2 and work rate) were small enough to separate the dynamic influence of each variable. Multivariate autoregressive analysis is an attractive method because the influence of one variable on another can be estimated and the cause and the result can be discriminated. However, as the basic principle of this method is to fit a linear model to time series data, there are some limitations in applying this method: (i) If certain innovations, i.e. residual remaining after the deletion of the effect of its past history, show strong correlations, their simultaneous inclusion in the model must be avoided. (ii) D a t a must include enough frequencies to determine the dynamic characteristics. 54

(iii) For a small magnitude forcing, the nonlinear system can be linearised about the operating point. The autoregressive model provides information on linearised system behaviour. However, it may only be good for the operating point studied and the magnitude of the forcing used. There is some evidence to suggest that the contributions of nonlinearities to the ventilatory response are small when the work intensities are moderate (BENNETT et al., 1981). However, it is well known that the C O 2 control system has nonlinear dynamics even though the steady-state fee - P a C O z response is linear (GRODINS et al., 1954). Moreover, the ventilatory response to combined hypercapnia and exercise has been reported to have a multiplicative component (CLARK et al., 1980; DUFFTN et al., 1980; HULSBOSCH et al., 1981) even in a steady-state, although other investigators have argued for a simple addition

Medical & Biological Engineering & Computing

January 1992

(BRADLEY e t al., 1980; KELLY et al., 1982). The failure of the superposition principle in our results clearly shows the nonlinear system dynamics. However, the present method does not discriminate the nonlinear characteristics of the controller from those of the plant. We have examined the responses to the following: (a) Random forcing of CO2 in the inhaled gas between 0 and 7 per cent. (b) Random forcing of work rate between 20 and 80 W. (c) The combination of (a) and (b). For a linear system, the operating point is defined by applying a constant stimulus equal to the mean of the random sequence. For the work rate sequence used in this study, that would be 50W. It is clear that the operating points for the three experimental conditions would all be different. Thus, the linearised models obtained would be characterisations of the system over three possibly different operating ranges. Therefore, the decrease in the relative contribution of the work rate during simultaneous C O 2 and exercise forcing may reflect the different operating ranges between the three experimental conditions. It is also possible that the actual magnitude of the C O 2 stimulus is different between rest and exercise. During exercise, the ventilation is higher, and thus the airway load of CO2 is also greater, resulting in a greater change in P a C O 2. Moreover, the alveolar-arterial Pco2 gradient increases with work load (DENISON et al., 1971). Although these effects are taken into account to some extent in the modelling, it means that the magnitude of the forcing signal has been changed. Thus, it will affect the estimation of AR coefficients and the partitioning of the linearised power spectrum. The results of the relative power contribution also showed that ventilation is largely unaffected by work rate and end-tidal CO2. The large contribution of the residual power may suggest that the variability of ventilation at this frequency range is composed of not only stochastic noises but also certain determininstic factors. BENNETT et al. (1981) showed that, to detect the rapid ventilatory ressonse, the system must be excited over a range of frequencies greater than 1.7 cycle min -1. As the present study uses a pseudorandom input with a 3 dB bandwidth of 1 cycle m i n - 1, our results do not excite a rapid response over 1 cycle m i n - t . However, that cannot explain the low coherency at 1 cycle rain -1 because the inputs we used have enough power to excite the system at 1 cycle m i n - i . Several explanations for the fluctuation of ventilation around this frequency can be put forward: (i) Contributions by the fluctuation of unobserved inputs (metabolic rate for C O 2 and 02 etc.). (ii) Variations independent of chemical feedbacks produced by central pattern generator. (iii) Behavioral, nonhomeostatic components. The decrease in the residual power when exercise is combined with CO2 inhalation might reflect the increased homeostatic demands relative to behavioural demands. In conclusion, we have shown that the dynamics of the response to C O 2 inhalation, exercise and their combination are nonlinear and that the combination of CO2 inhalation and exercise magnifies the nonlinear behaviour. Ventilation is largely unaffected by either work rate or end-tidal CO 2 at 1 cycle min-1. During simultaneous CO 2 and work rate forcing, ventilation tends to follow the change in the end-tidal CO2. However, our method has not specifically identified whether these nonlinearities represent characteristics of the controller or the plant. A better application of autoregressive modelling to obtaining

Medical & Biological Engineering & Computing

the controller characteristics might be a comparison of dynamic responses among subjects under the same operating point and the same magnitudes of input signals.

Acknowledgments--This work was supported in part by the ISM Co-operative Research Program (87-ISM-CRP-77).

References AKAIKE,H. (1967) Some problems in the application of the crossspectral method. In Spectral analysis of time series. HARRIS,B. (Ed.), John Wiley, New York, 81-107. AKAIKE, H. and NAKAGAWA, T. (1972) Statistical analysis and control of dynamic systems. Saiensu-sha, Tokyo, (In Japanese, with a computer program package TIMSAC written in FORTRAYIV with English comments.) BENCHITRIT,G. and BERTRAND,F. (1975) A short-term memory in the respiratory centres: statistical analysis. Respir. Physiol., 23, 147-158. BENNETT,F. M., REISCHL,P., GRODINS,F. S., YAMASHIRO,S. M. and FORDYCE, W. E. (1981) Dynamics of ventilatory response to exercise in humans. J. Appl. Physiol., 51, 194-203. BRADLEY, B. L., MESTAS,J., FORMAN,J. and UNGER, K. M. (1980) The effect on respiratory drive of a prolonged physical condition program. Am. Rev. Respir. Dis., 122, 741-746. CASABURI,R., WntPP, B. J., WASSERMAN,K., BEAVER,W. L. and KOYAL, S. N. (1977) Ventilatory and gas exchange dynamics in response to sinusoidal work. d. Appl. Physiol., 42, 300-311. CLARK, J. M., SINCLAIR,R. D. and LENOX,J. B. (1980) Chemical and nonchemical components of ventilation during hypercapnic exercise in man. J. Appl. Physiol., 48, 1065-1076. DAUBENSPECK,J. A. (1973) Frequency analysis of CO 2 regulation: afferent influences on tidal volume control. J. Appl. Physiol., 35, 662-672. DENISON, D., EDWARDS, R. H. T., JONES, G. and HOPE, H. (1971) Estimates of the CO 2 pressures in systemic arterial blood during rebreathing on exercise. Respir. Physiol., 11, 186-196. DUFFIN, J., BECHBACHE,R. R., GOODE, R. C. and CHUNG, S. A. (1980) The ventilatory response to carbon dioxide in hyperoxic exercise. Respir. Physiol., 40, 93-105. GRODINS, F. S., GRAY, J. S., SCHROEDER,R., NORINS, A. L. and JONES, R. W. Respiratory responses to CO 2 inhalation. (1954) A theoretical study of a non linear biological regulator. J. Appl. Physiol., 7, 283-308. HULSBOSCH,M. A. M., BINKHORST,R. A. and FOLGERING,H. T. (1981) Effects of positive and negative exercise on ventilatory CO 2 sensitivity. Eur. J. Appl. Physiol. Occup. Physiol., 47, 73-81. KELLY, M. A., OWENS, G. R. and FISHMAN, A. P. (1982) Hypercapnic ventilation during exercise: effects of exercise methods and inhalation techniques. Respir. Physiol., 50, 75-85. LAMARRA, N., WHIPP, B. J., WARD, S. A. and WASSERMAN,K. (1987) Effect of interbreath fluctuations on characterizing exercise gas exchange kinetics. J. Appl. Physiol., 62, 2003-2012. POON, C. S. and GREEN, J. G. (1985) Control of exercise hyperpnea during hypercapnia in humans. J. Appl. Physiol., 59, 792797. SOHAB, S. and YAMASnmO,S. M. (1980) Pseudorandom testing of ventilatory response to inspired carbon dioxide in man. J. Appl. Physiol., 49, 1000-1009. SUYAMA, A., SUNAGAWA, K., HAYASHIDA, K., SUG|MACHI, M., TODAKA, K., NOSE, Y. and NAKAMURA, M. (1988) Random exercise stress test in diagnosing effort angina. Circulation., 78, 825-830. SWANSON, G. D. and BELLVILLE, J. W. (1974) Hypoxichypercapnic interaction in human respiratory control. J. Appl. Physiol., 36, 480~87. WASSERMAN, K., WHIPP, B. J., KOYAL, S. N. and BEAVER,W. L. (1973) Anaerobic threshold and respiratory gas exchange during exercise. J. Appl. Physiol., 35, 236-243. WIGERTZ, O. (1977) Dynamics of ventilation and heart rate response to sinusoidal work load in man. J. Appl. Physiol., 29, 208-218.

January 1992

55

Appendix Estimation of the relative power contribution usin 9 a multivariate autoregressive model If the number of variables observed for the system concerned is K, we can obtain the K x N data matrix X(s), s = 1, 2 . . . . . N, during N sampling times. Their mean values are subtracted from the original data beforehand. For a system that is linear, time invariant and has a finite memory, the relationship between input xj(s) and output xi(s) can be written as K

xi(s) =

M

Z

~ Ai~(m) . xj(s -- m) + ui(s)

j=l,i~j

m=l

( i = 1,2 . . . . . K)

(1)

where xi(s ) is the output, xj(s); j = 1, 2. . . . . K, are the inputs of the system and ui(s) is the noise of x~(s) after deleting all influences from xj(s). Aii(m) is the impulse response function from the input xj to the output xi. If the variables include exogenous variables such as the work rate, regression coefficients to these variables are set to zero because they are not influenced by any other variable. The correlation between ui(s) and xj(s - m) is minimised when the noise is represented by its innovation, i.e. the residual remaining after the deletion of the effect of its past history. Whether there is feedback or not, this method is useful to whiten nonwhite or structured noise. Indeed, it is likely that the noise term is not white; its source remains unclear, although some investigators have proposed the involvement of an element of short-term memory (BENCHITRITand BERTRAND, 1975; LAMARRA et al., 1987). If we assume noncorrelation between innovations, the unbiased maximum likelihood estimates of Ai~(m ) can be calculated, i.e. impulse response functions and frequency response. M is determined as the order that minimises Akaike's information criterion (AIC) which is defined as follows: AIC = - 2 In (maximum likelihood) + 2(number of parameters) In this case, model order M can be determined objectively to preserve as much information as possible and to minimise the number of variables9 The system order M is not necessarily identical for all inputs. By assuming orthogonality between elements of the innovation, the relative power contribution of uj(s) to ui(s) can be computed, i.e. the power spectrum of x~(s) is decomposed into the sum of the contributions from xj(s) as follows: Aq(f), the frequency response function from u~(s) to ui(s), is given by M

Aii(f) = ~, Aij(m) exp (-i2ufm)

(2)

m=l

Under the assumption of orthogonality between elements of the innovation, Pu(f), the power spectral density function of each variable xi(s) is decomposed into the sum of the contributions from p(uiXf), the power spectral density function of uj(s): k

p.(f) = ~ Ibq(f) 12p(uj)(f)

(3)

j=l

where b~j(f) is the closed-loop frequency response function from xj(s) to xi(s ) and is expressed by

b.!:))

b i t ( f ) bl2(f).., blk(f)\ b2,.(f) b22.(f) 999

b.)S)..,

.

,iS)l b 1 -- alx(f)

__-

56

--a,2(f) ...

--a21(f) I

--

--akl(f)

--ak2(f)

aee(f )

--axk(f)\-' -

-

1 --

l:k(f) I a,,(S)l #

(4)

Therefore, the relative contribution of the variable xj to Pu(f) at the frequencyfis described by

qq(f) rij(f) = Pii(f)

(5)

where

qij(f) = ]bij(f)[2p(uj)(f)

(6)

Analysis of the correlation coefficients between the elements of innovation is useful for the selection of system variables to be used eventually. If certain innovations show strong correlations, considerable errors are introduced in the estimation of the relative power contribution and thus their simultaneous inclusion in the model must be avoided.

Authors' biographies Yoshitaka Oku was born in 1959 in Sumoto, Japan. He received his MD in 1983 from Kyoto University. He is currently a Research Associate of the Department of Medicine at Case Western Reserve University. His research interests are the control of breathing, especially rhythm generation, pattern formation and its modulation by chemical and neuromechanical inputs. Kazuo Chin was born in Osaka, Japan, in 1955, but he is Korean9 He received the MD degree in 1981 from the University of Kyoto, Kyoto, Japan. He is an instructor at the Department of Clinical Physiology, Chest Disease Research Institute, Kyoto University, Japan9 His current research interest is respiratory regulation and abnormal breathing during sleep. Michiaki Mishima was born in Osaka, Japan, in 1951. He received the MD degree in Internal Medicine in 1977 from the Kyoto Univeristy, Kyoto, Japan. He is an Assistant Professor at the Faculty of Clinical Physiology, Chest Disease Research Institute, Kyoto University, Japan. His current research interests include digital signal processing and modelling of the respiratory system. Motoharu Ohi was born in Saitama, Japan, in 1948. He received the MD degree in 1973 from the University of Kyoto, Kyoto, Japan9 He is an Associate Professor at the Department of Clinical Physiology, Chest Disease Research Institute, Kyoto University, Japan9 His clinical interests are regulation of respiration, abnormal breathing during sleep and exercise physiology. Kenshi Kuno was born in Kumamoto, Japan, in 1933. He received the MD degree in 1960 from Kyoto University, Kyoto, Japan. He is a Professor at the Department of Clinical Physiology, Chest Disease Research Institute, Kyoto University, Japan. His current research interest is respiratory insufficiency.

Yoshiyasu Tamura was born in Kobe, Japan, in 1952. He received the Ph.D. degree in Physics in 1980 from the Tokyo Institute of Technology9 He is an Associate Professor at the Statistical Data Analysis Center, the Institute of Statistical Mathematics. Since 1981 he has been engaged in time series analysis9 His current research interest is artificial intelligence, neutral networks, fuzzy relations and computer statistics9

Medical & Biological Engineering & Computing

January 1992