Med. & Biol. Eng. & Comput., 1978, 16, 262-268

'Estimation of frequencies of gastrointestinal electrical rhythms using autoregressive modelling D. A. Linkens

S.P.

Datardina

University of Sheffield, Dept. of Control Engineering, Mappin Street, Sheffield $1 3JD, England

method of determining the periodicities of gastrointestinal data using an autoregressive modelling technique is presented. The method has been applied both to simulated sinusoidal data and to gastrointestinal and colonic data. Unlike the fast Fourier transforms, the method presented is both automatic and capable of yielding frequency estimates on as little as 5 cycles of data corrupted with noise.

Abstraet--A

Keywords--Autoregressive mode], Data analysis, Gastrointestinal spectral analysis

1 Introduction

IT is well known that spontaneous electrical rhythms can be recorded from many parts of the mammalian gastrointestinal tract (DuTmE, 1974). The physiological mechanism responsible for the occurrence of this activity is basically unknown, and it has therefore not been possible to develop an empirical mathematical model of the type developed by Hodgkin and Huxley to model the behaviour of the squid axon. It has been found, moreover, that the use of extracellular electrodes and the presence of artefacts have meant that neither the actual shape nor the amplitudes of the electrical poteotial have contributed significantly in the modelling of this activity. However, the measurement of the frequency of the signals along the gastrointestinal tract has been possible, and this has indicated the presence of a frequency gradient along the tract. This has resulted in the hypothesis of a model comprising a series of bidirectionally coupled Van der Pol type oscillators for both the small and the large intestine (NELSEN and BECKER, 1968; ROBERTSON-DUNN and LINKENS, 1947; LINKENS et al., 1976) and a 2dimensional array of such oscillators for the stomach (SARNA e t al., 1972). The gastrointestinal frequencies have normally been estimated by the application of fast Fourier transforms (f.f.t.) on stretches of data (LINKENS and CANNELL, 1974) but this has a number of disadvantages. To obtain sufficient frequency resolution the f.f.t, requires fairly long stretches of data covering many complete cycles. Thus a stretch of periodic gastric data with a frequency of 0' 2~3.3 Hz has to be about eight minutes long to yield a reasonably accurate estimation of frequency, and Received 20th August 1976

this is the maximum length for which the data may be normally taken to be stationary. Also because of the noisy nature of the signals it is sometimes difficult to determine either visually or automatically which are the significant peaks in an f.f.t, plot. It is therefore of value to develop a method of measurement that is capable of direct and accurate estimation of frequency on shorter stretches of data. In the in vivo experimental situation a minicomputer provided with an analogue/digital convertor linked to the experimental subject would enable the investigator to build up a comprehensive picture of the changes in frequency with time in various parts of the alimentary canal. Stationary e.e.g, signals have been modelled by a linear difference equation with constant coefficients. The parameters of the equation have been estimated by autoregressive moving-average techniques (ZETTERBURG, 1969). More recently, by employing only the autoregressive technique, these signals have been modelled by a linear filter having constant coefficients in the denominator only. The parameters of the filter are held to be appropriate if the residuals from the estimation are approximately white noise. The linear filter has been used (LoPEs DA SILVA et al., 1975) to detect spike and wave complexes typical of epileptic patients. This has been done by locating the regions in the e.e.g, where the output from the filter is not white noise and hence indicates that the signal is nonstationary. The filter has also been used to determine the average power and ecntre frequency of the =-activity (PFURTSCHELLER, 1976). In this communication the gastro-intestinal data is considered to be stationary over short stretches of time and auto-regressive modelling technique is applied to detect the periodicities in the data.

0140-0118/78/0618-0262 $1 950/0 9 IFMBE: 1978

262

Medical & Biological Engineering & Computing

May 1978

2 Theory

estimated parameter vector ~x...~p is given by

The method of autoregressive modelling applied to e.e.g, data has been described in detail by (DIJK, 1973) and (LOPES DA SILVAet al., 1975). The application of this method to gastrointestinal signals follows a similar course so that only a brief summary of the method will be presented here. A sampled gastrointestinal signal at any particular instant of time can be represented by a weighted sum of its own past together with an error (noise) term as follows:

y(k) = - ~ 1 y ( k - 1 ) - ~ 2 y ( k - 2 ) . . . . . . - ~,, y ( k - p ) + ~ ( k )

(1)

[il i, ...... =

,

:

7

:

~ r(p)_l

Lr(p - 1) ...... r(0)_]

(2)

where r(o).., r(p) are the estimated correlations of the stationary time series given by the first p lags of the correlation functions of the gastrointestinal signal. Inverting the matrix in eqn. 2 is not a major problem as the order of the matrix is only p x p. To determine if the residual is white noise an estimation of its variance is required. This can be shown to be t~2= r(0)+ ~ ~ir(i) i=1

Eqn. 1 represents a pth order auto-regressive (a.r.) model where ~(k) is the noise term with variance a 2, ctl ... Ctpare the parameters of the model, y ( k - 1) = y ( ( n - 1 ) T) for n = 1...N is the value of the sampled gastrointestinal signal at the ( k - 1 ) t h instant, N is the total number of points and T is the sampling period. The parameters ~q...~tp are estimated by minimising the sum of the squares of the residuals at each instant in time by a least-squares procedure. The 1"596E0

II

,7

-1596E0

Introducing the shift operator z so that y ( k - 2 ) = z-1 y(k) and substituting into eqn. 1 leads to y(k)(1 +~q z -a + ... +r

-~) = ~(k)

(4)

Eqn. 4 represents a linear filter in transform variable z. The roots of the polynomial in eqn. 4 correspond to the poles of the filter. A pair of poles found to be on the unit circle in the z-plane implies that the data is periodic. Multiple pairs of poles around the unit circle indicate either several periodicities or harmonic components, and knowing the sampling period the frequencies of these components can be easily computed. A pair of poles is considered to be significant if the distance of each pole from the origin is greater than 0.95 and less than 1.0 (an empirical number based on results derived from typical stretches of gastrointestinal data.) Fig. lb shows a typical result with two pairs of poles on the unit circle indicating the presence of a fundamental and the 2nd harmonic. The values of frequencies corresponding to these poles can be printed out directly. 3 D a t a collection from digestive tract

Fig. 1A Computer-generatednoise-corruptedsine-wave plus second harmonic. The signal has been scaled horizontally and vertically

3.1 Experimental methods The in vivo human data analysed in this work was obtained by means of an intraluminal suction

Zeros

Poles Real 9' 378E-- 1 9"378E-- 1 - - 8" 0 2 5 E - - 1 --8"025E--1 7" 8 8 4 E - - 1 7" 8 8 4 E - - 1 --2" 813E--1 --2"813E--1

Imaginary 2" 9 6 3 E - - 1 - - 2 . 963E - - 1 --3" 702E-- 1 3"702E--1 - - 5" 6 3 2 E - - 1 5" 6 3 2 E - - 1 --8"075E--1 8"075E--1

Real O" O00E-- 7 O" O 0 0 E - - 7 O" O00E - - 7 O'O00E--7 O" O00E-- 7 O- O00E - - 7 O'O00E--7

Imaginary O" O00E-- 7 O" O 0 0 E - - 7 O" O00E - - 7 O-O00E--7 O" O00E - - 7 O" O 0 0 E - - 7 O'O00E--7

Fig. 1B Z-plane plot of the factorised polynomial of an 8th-order autoregressive model fitted to data in Fig. IA. The poles in the 1st and 4th quadrant are estimates of the fundamental and the second harmonic

Medical & Biological Engineering & Computing

May 1978

263

electrode mounted in the side wall of a tube 5 mm in diameter. The stainless-steel electrode was 0.25 m m in diameter and surrounded by a rubber cap which enabled the electrode to be sucked onto the inner wall of the bowel using a negative pressure of 20 cm Hg, This gave a mucosal penetration depth of about 2 mm. F o r the canine recordings serosal electrodes were used (WATERFALLet al., 1972) with the connections being taken through a surgical drain. The recordings were monopolar and were all amplified (frequency response 0.02 Hz to 1 kHz) before being recorded onto 4-track analogue magnetic tape. Further details of the recording methods can be found in TAYLORet aL, 1974. 3.2 Implementation The algorithm is implemented on a G E Compac 4020 using the Sheffield package for the analysis and identification of data (BATEY et al., 1975). The package incorporates about 30 options among which are the following facilities relevant to this paper: (a) generation of frequency a n d / o r amplitude modulated sine waves with noise (b) digital filtering of data (c) fitting an autoregressive model to the data (maximum order of 20) (d) computing the roots of polynomial (maximum order of 20) (e) computing and displaying the autocorrelation function of the residuals 4 Results

The following stretches of simulated and recorded data were analysed for periodicity: (a) simulated sine wave and noise (b) simulated sine wave and harmonic and noise (c) recordings from canine duodenum and stomach (d) filtered recordings from human stomach (e) recordings from human colon. F o r the case of a noise-contaminated signal comprising a sine wave plus a second harmonic component some typical results photographed from the computer screen are shown in Fig. 1. The data being analysed is scaled vertically and horizontally before being displayed as in Fig. la. After the autoregressive model coefficients have been found the resulting polynominal is factorised and the real and imaginary parts of the poles displayed as shown in Fig. lb. These poles are also plotted on a z-plane where the relevant spectral components in the signal can be easily recognised by their location in the first and fourth quadrant and their proximity to the unit circle. The fundamental and second harmonic frequencies from the poles of Fig. lb are 48.70 and 9 8 . 7 2 H z while the equivalent f.f.t, on 512 points of data gave 47.1 and 96.1 Hz. The frequency 264

discrimination, however, in this latter case was + 1" 95 Hz. In the following results obtained from canine and human digestive recordings the sampling frequency was set at 2 Hz. This has been found to be a convenient and useful value for most digestive signals although it is not optimum in terms of frequency discrimination using f.f.t. The frequency discrimination for a 512 point f.f.t, for all the following results was 0.008 Hz. A 100 point stretch of canine duodenal data is shown in Fig. 2a and the transform for 512 points gave clear spectral components at the fundamental and second harmonic frequency. Fig. 2b shows how significant components are calibrated using the screen cross wires and shows the low-frequency artefacts which often occur in digestive electrical recordings. F o r this data a 6th-order a.r. model was sufficient to give a good estimation as seen in the pole-zero plot of Fig. 2c

2.339E2

11

Li

i~ Fig. 2A 100 points of canine duodenal data sampled at a frequency of 2 Hz

fast fourier transform of channe| one ,2-210E3

[ s q u a r e of d.c.=5-686E3

B

100

200

-2210E3 A = 2. 941 E-- 1 B = 5.921E--1

cycles cycles

A M P = 2" 035E3 AMP = 7"447E2

Fig. 2B Spectral analysis of 512 points of canine duodenal data. The frequency discrimination is 0.008 Hz

M e d i c a l & Biological Engineering & C o m p u t i n g

M a y 1978

and verified by the residuals of Fig. 2d which were all within the + 1 standard deviation (s.d.) lines. In contrast to canine duodenal signals which have relatively small harmonic content the canine gastric recordings are pulse-like in nature. This can be seen in Fig. 3a where the rhythmic sharp spikes were due to the electrical basic rhythms and the large pulses following these were related to the

mechanical contractions. In attempting to fit an autoregressive model to this data it was found necessary to use a 12th-order model which gave a good estimate of the fundamental frequency as 0.0956Hz (0.0902Hz by a 512 point f.f.t.) and gave residuals outside the _+1 s.d. lines as shown in Fig. 3b. The human stomach also gives electrical signals

Poles Real 3.011E--1 --8.188E--1 --8"188E--1 9.835E--1 --2.518E--1 --2.518E--1 5.895E--1 5.895E--1

Zeros Imaginary 0.000E--7 3.795E--1 --3-795E--1 0.000E--7 --9-363E--1 9.363E--1 --7.884E--1 7-884E--1

Real 0.000E--7 O-000E--7 0.000E--7 0.000E--7 0.O00E--7 0.000E--7 O.O00E--7

Imaginary 0.000E--7 0.000E--7 0-O00E--7 0.O00E--7 0.000E--7 0.O00E--7 0.000E--7

Fig. 2C Z-plane plot of the factorised polynomial of a 6th-order a.r. model. The estimated fundamental and the second harmonic were O. 029 and O. 058 Hz

-1.000E0

, O 0 0 E O

*ls.d

.

I

.

.

.

~.... , . J ~ , A , , . . . . . . . . . A, ,.6'~,, A,, J~..,,, / v -ls.d[--

" l~v ---

V~o'-

~vu

~o ~ ,

,ls.d. ,

. . . . .

,f?~.,,, ~Q-~'-~..,, .....t~.. ,z~. ,AM,, ,

9

I-I'O00EO Fig. 2DAutocorrelated function of 50 lags on residuals of the canine duodenal data. As autocorrelation function is well within the standard deviation line, the residuals can regarded as white noise

the the +_1 be

Fig. 3B Autocorrelation function of 50 lags on the residuals of the canine gastric data. The residuals are from a 12th-order model

,7.739E2 -2'043E2

'1'1

.

,t~. ,=

-?.739E2

j j

200

L~J

-2"043E2 Fig. 3,4 200 points of canine duodenal data. Sharp pulses indicate electrical activity and the large after potentials indicate mechanical contraction

Medical & Biological Engineering & Computing

Fig. 4A 200 points of human gastric data showing square-like rhythms

May 1978

265

with high harmonic content as seen in the 200-point display of Fig. 4a. To fit this type of data it was found to be desirable to lowpass filter the data

several components spread around a centre frequency as seen in Fig. 5b. A 15th-order a.r. model fitted this data well and gave two spectral components near

Poles

/

Real 9-623E--1 2"484E--1 2"484E--1 4-470E--1 4"470E--1 --7"984E--1 --7"984E--1 9-302E--1 9-302E--1 7.992E--1 7'992E--1 --1-994E--1 --1"994E--1 --6"088E--1 --6-088E--1

Zeros Imaginary 0-000E--7 --7"823E--1 7"823E--1 7"430E--1 --7'430E--1 1"328E--1 --1"328E--1 --2"643E--1 2"643E--1 --5"051E--1 5-051E--1 8"161E--1 --8"161E--1 --6"128E--1 6-128E--1

Real O-000E--7 0-000E--7 0"000E--7 O-000E--7 0'000E--7 0'000E--7 0"000E--7 0-000E--7 0"000E--7 0-000E--7 0"000E--7 0-000E--7 0.000E--7 0.000E--7

Imaginary O-000E--7 0.000E--7 0"000E--7 0.000E --7 0.000E--7 0-000E--7 0"000E--7 0-O00E--7 0.000E--7 0.000E--7 0"000E--7 0.000E--7 0"000E--7 0.000E--7

Fig. 4B Z-plane plot of a 15th- order autoregressive model fitted to human gastric data that has been Io wpass filtered to remove components above the second harmonic

before modelling to remove components above the second harmonic. The pole-zero plot for a 15thorder model is shown in Fig. 4b where the fundamental- and second-harmonic components have been estimated at 0.0578 and 0.1161 Hz (0.055 and 0.11 Hz, respectively, by a 512-point f.f.t.). Even with this high-order model the residuals were slightly outside the + 1 standard-deviation lines in a manner that suggested that the fundamental had not been estimated perfectly. The fundamental, second and fourth harmonic can however be estimated by a 20th-order autoregressive model. The third harmonic was not estimated presumably because the power spectral component corresponding to the fourth harmonic has greater magnitude than that of the third harmonic. The last results shown are for human colonic signals which show much more variability in frequency. The f.f.t, of a 512-point stretch has

the unit circle as shown in Fig. 5c. The fundamental and second harmonic estimates were 0-0881 and 0.178 Hz, respectively, while the residuals lie well within the + 1 s.d. lines. Besides the stretches of data mentioned in (a) to (e) above, the method was also applied to the following simulated data with noise corruption: (i) frequency- a n d / o r amplitude-modulated wave forms (ii) composite waveforms not harmonically related. In eqn. 1 the fundamental could be identified with an 8th-order model but the side bands, clearly shown by an f.f.t., were not identified. In eqn. 2 2"364E2

-3.746E2 100

200

-2.364E2 = 5.882E--2 = 8.627E--2 = 1-019E--1 -3'746E2 Fig. 5A 250 p o i n ~ of human colonic data

266

cycles cycles cycles

AMP = 1.620E2 AMP = 2"364E2 AMP = 1 -775E2

Fig. 5B F.F. T on a 512-point stretch of colonic data. The frequency discrimination is O. 008 Hz

Medical & Biological Engineering & Computing

May 1978

both the periodicities of a composite waveform were detected with a 14th-order model. A number of heuristic models have been proposed to determine the optimal order of the a.r. models

University of Sheffield, Department of Control Engineering) for his help in the computer implementation of the algorithm, andProf. Duthie and his team at the Department of Surgery for furnishing the gastrointestinal data

Poles Real 3.357E 3. 357E 9. 035E 9. 035E

-----

1 1 1 1

Zeros Imaginary 6. 994E --- 6- 994E --- 3" 453E -3- 453E --

-- 6.938E-- 1 -- 8. 736E --- 8. 736 E ---8,027E--1 -- 8. 027 E -9. 505E -9- 505E-7.166E -7.166E ---4- 794E ---4.794E--1

2 2 1 1 1 1 1 1

1 1 1 1

O- O00E - - 7 -- 7. 043E -7- 043 E ---3.443E--1 3.443 E -- - 1 97 4 7 E - 1 97 4 7 E - -- 5. 972E-5. 972E --- 6' 999E -6.999E--1

1 1 1 1 1 1 1 1

7 7 7 7

Imaginary O. O 0 0 E - - 7 O, O 0 0 E - - 7 O. O 0 0 E - - 7 O, O 0 0 E - - 7

O. O 0 0 E - - 7

O. O 0 0 E - - 7

O- O00E - - 7 O, O 0 0 E - - 7 O-O00E--7 O. 0 0 0 E - - 7 O. O 0 0 E - - 7 O. O00E - - 7 O. O00E - - 7 O. O00E - - 7 O. O00E - - 7

O. O 0 0 E - - 7 O, O 0 0 E - - 7 O.O00E--7 O, O 0 0 E - - 7 0 9O 0 0 E - - 7 O. O 0 0 E - - 7 O. O 0 0 E - - 7 O. O 0 0 E - - 7 O. O 0 0 E - - 7

O. O. OO.

Real O00E-O00E-O00E-O00E--

Fig, 5C Z-plane plot of the 15th-order autoregressive model fitted to colonic data. The estimates of fundamental and second harmonic are 0-088 and O-178Hz

but these have not been found to be satisfactory (LOPES DA SILVA et aL, 1975). The optimal model orders quoted above were determined by tlial and error. Modelling a stretch of data with a higherorder model than necessary did not necessarily lead to greater accuracy. The model order required to give a good estimation could be reduced by filtering the data before estimation. The number of complete cycles necessary for analysis could be as little as two or three for simple data while about five cycles were required for more noisy data and harmonics content. 5 Conclusion

it is apparent from Figs. 1-5 that the autoregressive techniques applied to the determination of periodicities in gastrointestinal data is effective. The fundamental frequency can be fairly accurately determined even with a short stretch of data. To estimate the harmonics, longer stretches of data and higher-order models are necessary. The order of the models varies with the type of data being analysed. However, the order can be reduced by bandpass filtering the data so that higher harmonic components and very-low-frequency variations are removed. The method is now being transferred to a Data General Nova 16-bit machine with a 16k bit core, a back-up disc and Tektronix v.d.u, facility. Acknowledgments--The authors thank Dr. S. Billings, Medical

& Biological Engineering & Computing

on which this work has been based. S. Datardina acknowledges support of the UK Science Research Council in this research.

References BATEY, D. J., STERLING, M. J. H., ANTCLIFFE,D. J. and BILLINGS, S. A. (1975) 'The design and implementation of an interactive data analysis package for a process computer', Comput. Aid. Des., 7, 265-269. Box, G. E. P., and JENKINS, G. M. Time series analysis, ]brecasting and control, chap. 3, Holden Day, 1969. D1JK, A. (1973) Estimation of the dynamic properties of the e.e.g, by means of linear models. Report, Laboratory of Technical Physics, Technical University, Delft, (in Dutch). DtJTHIE, H. L. (1974) 'Electrical activity of gastrointestinal smooth muscle', Gut, 15, 669. LINKENS, D. A., and CANNELL,A. E. (1974) 'Interactive graphics analysis of gastrointestinal electrical signals', IEEE Trans., BME-20, 335-339. LINKENS, D. A., TAYLOR,I., and DUTHIE, H. L. (1976) 'Mathematical modelling of the colorectal myoelectrical activity in humans', Ibid., BME--22, 101110. LOPESDE SILVA, F. H., DIJK, A., and SMITS,A., (1975) 'Detection of nonstationarities in EEG's using the auto-regressive model--an application to EEG's of epilectics'. Brain Research Department, Institute of Medical Physics TNO, Da Costakade 45, Utrecht. NELSON, T. S. and BECKER,J. S. (1968) 'Simulation of the electrical and mechanical gradient of the small intestine', Am. J. Physiol. 214, 749-757. May 1978

267

]aFURTSCHELLER, G. (1976) 'Spectral parameters estimation of ~z-activity by means of the auto-regressive model', Biomed. Comput. Tech. Foundation Issue. ROBERTSON-DONN, B. and LINKENS, D. A. (I974) 'A mathematical model of the slow-wave electrical activity of the human small intestine,' Med. & Biol. Eng. 12, 750. SARNA, S. K., DANIEL,E. E. and KINGMA,Y. J. (1972) 'Simulation of the electrical control activity of the stomach by an array of relaxation oscillators. Am. J. Dig. Dis. 17, 299.

268

TAYLOR,I., ])UTHIE,H. L., SMALLWOOD,R. and LINKENS, D, A. 1974, 'The effect of stimulation on the myoelectrical activity of the rectosigmoid in man'. Gut 15, 599-607. WATERFALL,W. E., BROWN, B. I-[., DUTHIE,H. L. and WHITTAKER,G. E. (1972) 'The effect of humeral agents, on the myoelectrical activity of the terminal ileaum'. ibid. 13, 528-534. ZEYrERBURG,L. H. (1969) Estimation of parameters for a linear difference equation with application for e.e.g. analysis. Math./Biosei., 5, 227-275.

Medical & Biological Engineering & Computing

May 1978