IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MAY 1976

262 1.0

,....,...

rz .? .0

(r&l

--

...-

O.? r

700

.4..

J*

-a

0-?1

0

.DI

-do

-03

.04

.*or

.09

.011

.08

oq

I

lo

fI

f2

.1F-ALS f-I RAf FALSE POSITIVE RATE

Fig. 9. Curves of detection rate versus false positive rate for ten control groups of 42 normal and 33 infarcted subjects using reference waveforms derived from ten randomly chosen test subsets of 42 normal and 33 infarcted subjects. The control group in each case is the set of subjects remaining in the total pool of 84 normal and 66 infarcted subjects after choice of the test set.

facilities, and Dr. G. Weiss of the National Institutes of Health for his advice and encouragement. REFERENCES 1 C. V. Nelson, "Human thorax potentials," Ann. N. Y. Acad. Sci., vol. 65, pp. 1014-1050, 1957. [21 D. B. Geselowitz, "Multipole representation for an equivalent cardiac generator," Proc. IRE, vol. 48, pp. 75-79, 1960.

[3] B. Taccardi, "Distribution of heart potentials on the thoracic surface of normal human subjects," Circ. Res., vol. 12, pp. 341-352, 1963.

[41 L. G. Noran, N. C. Flowers, and D. A. Brody, "Body surface

[5] [61

171 [8] [9] [101

[111 [121 [13] [14]

potential distributions: Comparison of naturally and artificially produced signals as analyzed by digital computer," Circ. Res., vol. 13, pp. 373-387, 1963. J. P. Boineau, M. S. Spach, T. C. Pilkington, and R. C. Barr, "Relationship between body surface potential and ventricular excitation in the dog," Circ. Res., vol. 19, pp. 489-496, 1966. R. H. Selvester, J. C. Solomon, and T. L. Gillespie, "Digital computer model of a total body electrocardiographic surface map: An adult male-torso simulation with lungs," Circ., vol. 38, pp. 684-690, 1968. R. A. Helm and T. C. Chou, "Computation of a variable location dipole representation from body surface leads," Amer. Heart J., vol. 77, pp. 363-366, 1969. J. H. Holt, A. C. L. Barnard, M. S. Lynn, and P. Svendsen, "The study of the human heart as a multiple dipole source I: Normal male subjects," Circ., vol. 40, pp. 687-696, 1969(a). J. S. Abildskov et aL, "Distribution of body surface potentials with experimentally-induced multiple cardiac generators," Adv. Cardiol., vol. 10, pp. 69-76, 1974. E. J. Fischmann, M. Haghighat, and P. P. Mehrotra, "Surface potential profiles in myocardial infarction," in Advances in Cardiology, S. Rush and E. Lepeschkin, Eds. Basel, Switzerland: Karger, 1974(a), vol. 10, pp. 275-282. F. Kornreisch and D. Brismee, II, "Diagnosis of left ventricular hypertrophy and myocardial infarction from total surface waveform information," Carc., vol. 49, pp. 996-1004, 1973. J. Herson, M. L. McNeel, G. H. Weiss, and E. J. Fischmann, unpublished observations. T. Kadota, "Optimum reception of binary sure and Gaussian signals," Bell Syst. Tech. J., vol. 44, no. 8, pp. 1621-1648, 1965. H. L. Van Trees, Detection, Estimation and Modulation Theory. New York: Wiley, 1968, p. 274.

A Flexible Renewal Process Simulator for Neural Spike Trains BARRY F. BELANGER, WILLIAM J. WILLIAMS, SENIOR MEMBER, IEEE, AND TOM C. T. YIN, MEMBER, IEEE Abstract-Certain neuronal spike trains may be viewed as stochastic, nonhomogeneous point process. Neuronal information may be encoded in the time-varying mean rate of the spike train in some cases. For this purpose the simulator can be modulated by arbitrary timevarying voltage waveforms.

INTRODUCTION The term renewal process is applied to a broad category of random processes wherein the time intervals between discrete events are independent, identically distributed random vanables. Such processes are observed in a wide range of physical and biological settings. The Poisson process is one well-known renewal process. In many cases the spike discharge of a neuron can be characterized as a renewal process if each action potential is taken to be the occurrence of an event. This is particularly true for neurons in the central nervous system. Recent reports sum-

marized research findings related to statistical analysis of neuronal spike data [1], [2]. Earlier hopes that such spike train data could provide structural information about the neural network involved have proven to be overly optimistic.

Manuscript received April 22, 1974; revised May 29, 1975. This research was supported in part by USPHS Grant NS 08470. B. F. Belanger and W. J. Williams 'are with Bioelectrical Sciences Laboratory, Department of Electrical and Computer Engineering and the Bioengineering Program, University of Michigan, Ann Arbor, MI 48109. T. C. T. Yin was with the Bioelectrical Sciences Laboratory, Department of Electrical and Computer Engineering and the Bioengineering Program, University of Michigan, Ann Arbor, MI 48109. He is now with the Department of Physiology, The Johns Hopkins University, Baltimore, MD 21205.

COMMUNICATIONS

The analysis and modeling of the basic process are still of interest, however. Much of the experimental and analytical work has focused on homogeneous renewal processes. The statistical parameters of a homogeneous process are fixed in time and do not change from interval to interval. More recently, several investigators have focused their attention on neuronal spike trains with modulated mean rates [3], [4]. A sinusoidally modulated neuronal Poisson process has been observed and analyzed [ 3 ]. The difficulty in obtaining mathematical expressions for the time autocorrelation and spectral density functions of renewal processes is considerable. Attempts to analyze more complex situations, i.e., modulation of generalized renewal processes, are quickly bogged down by mathematical complexities. We are studying the properties of neurons in the central nervous system that respond to position and movement of the limbs [111. Some of these neurons exhibit interspike interval statistics that are well described by exponential or gamma distributions. The firing rate of the neuron can be modulated periodically by applying sinusoidal motion to the limb. In order to test hypotheses on the nature of the processing of these modulated neuronal spike trains at various levels in the nervous system, we have devised a system for the simulation of modulated neuronal spike trains with interspike interval distributions having exponential and gamma properties. In addition, we have found the simulator to be of value in testing new computational techniques. It is often valuable to have a simulated spike train with known properties for use in evaluating a system intended for processing of neuronal spike data. Our simulator is not a model or neural analog. It represents an attempt to precisely simulate the mathematical description rather than model a biological entity that approximates the condition. Therefore, no physiological significance should be attached to the system that produces the modulated spike train. We can regard the spike train produced by the device as an idealization of a neuronal spike train. The device has been of considerable value in our work. GENERATION OF PSEUDORANDOM RENEWAL PROCESSES A pseudorandom sequence generator is the basis for the simulated renewal process. This approach was suggested by the all-or-none digital nature of the point process itself and by the desire to have flexible interarrival time probability density functions (PDF's) and a voltage-controlled mean rate. In our case the basic generator is a 15-stage shift register with linear, modulo-two feedback (see Fig. 1). The choice of feedback and number of stages is discussed by Korn [5] and Elspas [6], along with the statistical properties resulting from such arrangements. At any particular instant each flip-flop is in either a zero or one state, its next state being determined by the present state of its input. The clock pulses trigger the state changes through the shift register. This particular arrangement shown in the upper portion of Fig. 1 has three important properties. 1) The sequence is periodic (i.e., pseudorandom) of length P = 215 - 1 clock pulses. 2) The occurrence of a zero or one at any stage is equally likely. 3) The time intervals between one's appearing at any stage are independent within the period "P" of the shift register. The output of flip-flop 15 may be described as a random telegraph wave and has the three properties described above. By performing the logical AND operation with this wave and

263

Fig. 1. Block diagram of the electronic simulator. Insets A, B, C, D, E show typical voltage signals observed on the indicated lines. See text for details.

the clock pulse train, one obtains a pseudorandom pulse train [Fig. l(c)] with the three pseudorandomness properties. The statistics of this pulse train can be derived from the properties of the shift register. The clock pulses can be considered to be independent Bernoulli trials during which a pulse may occur with probability 1/2. The statistics are binomial and the counting function distribution is given by Prob[k successes in n trials] = (W)pkqflk where the (1) probability of a success (one) is "p" and the probability of a failure (zero) is "q." Here, p = q = 1/2. From this one can obtain the interarrival time distribution. Let

PISI(k) = Prob[(k - 1) failures in (k - 1) trials and a success in the kth trial]

PISI(k) = qk-l p PISI(k) = (1/2)k

(2)

(3) (4)

which is a geometric series in "k" (see Fig. 2). The Poisson process may be obtained from a series of Bernoulli trials with binomial statistics by letting "n," the number of trials, approach infinity while "p," the probability of success, approaches zero under the constraint that their product X remains finite and nonzero (7). Under this condition the binomial counting function converges to the Poisson

Prob[k successes in time (tI, t2 )] e

x(t2-tl)[X(t2

-

tI

)I k

t2 >tl.

(5)

The corresponding interarrival time PDF p (t) may be obtained as follows. Let P[t, At] = Prob[no success in (0, t) and one success in (t, t + At)] = Prob[no success in (0, t)] Prob[one success in (t, t + At)] from the independence property, -=t . e-XAt XAt where 0 < At 0) by ANDing the pulses with the output of a high-frequency multivibrator. In turn, the clock rate must be increased (n -> oo) in order to keep the mean rate of the final pulse train "E" at a finite nonzero value. The trial is seen to be independent through three considerations. 1) The multivibrator frequency is high compared to the clock rate and no harmonic relation exists between them. 2) The pulse width at "C" is small compared to the up-time of the multivibrator. 3) The probability of success is determined by the ratio of the up-time to the total period of the multivibrator. This operation cannot, of course, generate a true Poisson process, but the degree of convergence (p 0.07) is sufficient to yield an ISIH which closely approximates the continuous exponential [Fig. 2(b)], as compared to the histogram obtained from the system without attempted convergence -

[Fig. 2(a)J.

order gamma density with a logical operation which emits one pulse for every two pulses in the simulated Poisson process. Higher order gamma functions can be obtained easily by the addition of more flip-flops in the logic. The convergence imposed by the Bernoulli sampling process does not alter the fact that the interarrival times are independent, and in addition, the sequence is no longer periodic or repeatable. Fig. 3 shows serial correlations corresponding to pulse trains with discrete, exponential and gamma PDF's. The small correlation coefficients bear out the assumptions made above. MODULATION OF THE MEAN RATE The mean rate of the pulse generator X is dependent on the clock rate and the success probability at any clock pulse

An nth-order gamma function can be viewed as an n-fold convolution of an exponential function [8]. This can be achieved by summing n independent exponentially distributed random variables. If we consider the interpulse arrival times to be random variables, we know that they are independent and have identical exponential PDF's. Thus, to form an nth-order gamma interarrival time PDF one must sum n consecutive interevent intervals in a Poisson process. This amounts to summing n pulses in our case. We have formed the second-

(8)

where X = the mean rate of the random process. With the probability of success fixed, the mean rate is directly proportional to the clock frequency. If the clock used is a voltage-controlled pulse oscillator that can be frequencymodulated by arbitrary time varying voltage inputs, then the mean rate of the pseudorandom process is in turn modulated by the input signal. If nf(t) =no + nm * fm(t) (9) where n (t)

the instantaneous clock frequency the modulation parameter fm (t) the modulating signal then X(t)=p *n(t)=p -no pnm *fm(t) which may be written X(t)=X0 +Xmfm (t) where { xo = p no nm

Xm =p

(10)

(1 1)

nm.

It may be noted that the clock modulation depth nmInO remains unchanged in the renewal process, since

Xm

P nm nm (12) p no no This arrangement offers precise control over statistical and modulation parameters. In order to deal with the general case of modulated renewal processes it is necessary to relax the condition that event occurrences are identically distributed probabilistically. Instead, we incorporate the modulation into a parametric change in the interarrival time PDF. For instance, the Poisson process has an exponential interarrival time PDF with parameter X, the mean rate. A nonhomogeneous Poisson process is one in which X is allowed to be time varying or, equivalently, modulated. These nonhomogeneous renewal processes are thus a generalization of homogeneous renewal processes. Instead of expressing the properties of a nonhomogeneous renewal process in terms of the interarrival PDF it is often convenient to deal instead with the autocorrelation of the process [3] or the power spectral density (PSD) of the process. The autocorrelation and PSD are Fourier transform pairs. For periodically modulated renewal processes the PSD form is particularly attractive. The frequency components of the modulating waveform appear in the PSD at discrete frequencies identical to the frequencies present in the modulating waveform (see Fig. 4).

X0

COMMUNICATIONS

265

-f ui

OF

v."

4V\y/

la m

.

^ A

0.O

f.oW

-

so

cr__,

el..

a...s a....

.0

(a)

2.0

4.0

6.0

8.0

FREQUENCY (HZ)

10.0

~.

0 a..***

12.0

114.0

Fig. 4. PSD for simulated and neural trains. PSD obtained from a gracile nucleus unit responding to sinusoidal knee motion (circles) and the simulator driven by a sinusoidal waveform at the same frequency (solid line). The statistical and modulation parameters of the simulator were matched to those of the gracile cell unit. The correspondence is close and the modulating signal at 4.0 Hz is clearly evident in both PSD's. The discrepancy at low frequencies may be due to effects of respiratory rhythm.

Li

r

V

Li.

vV

/V C:~ " a." vOR O

'4 A

A

-

F

oo

'I.

(b)

U

6-4

LA.

8 1-

1

,

.

4

00

OIREO

v

V

,00,

A

51 8 (c)

Fig. 3. Serial correlation resultss. Serial correlations obtained from simulated pulse trains with ([a) discrete, (b) exponential, and (c) second-order gamma ISIH's. The small correlation coefficients are a necessary condition for inde-pendent intervals, as explained in the text.

APPL:ICATIONS The specific application of the simulator in our laboratory is in a study of the coding andI transmission of joint position information. At the level of t;he nucleus gracilis (the secondary

neuron location in the relay chain) neuronal firing patterns are sometimes observed with characteristics similar to modulated renewal processes. A description of the experimental approach and results is available [ 11 ] . When the knee joint of the experimental animal is articulated with a joint angle waveform that is sinusoidal in time, the neuronal firing pattern at the gracile level sometimes exhibits PSD characteristics quite similar to the PSD characteristics obtained with the simulator. Fig. 4 compares the PSD's obtained experimentally and with the simulator. Here the average firing rate and the depth of modulation of the simulator were adjusted to match the experimental data. The PSD's then obtained compare quite favorably. Modulated renewal process firing patterns can be quite deceptive to the ear, especially with low modulation depths. Often one is not able to detect the sinusoidal component in the erratic firing pattern. Since one often depends on sound patterns to guide the search for single neuronal units it is difficult to isolate neurons with renewal type firing patterns responding to stimuli. In spite of the erratic sound of the firing pattern, the autocorrelation and PSD are often quite consistent and the component at the modulating frequency is clearly evident. The simulator is useful in several ways as an aid to interpretation of experimental data. Some neuronal firing patterns produce gamma type ISI's rather than exponential ISI's. Modulated gamma processes are easily produced by the simulator for comparison with experimental results. The simulator can also be driven by triangular or square modulating waveforms for comparison with experimental results obtained with such inputs. Time delay and refractory period are easily incorporated into the simulator in a controlled fashion. To account for all of these conditions analytically is a difficult undertaking. In addition to simulation of modulated neuronal firing patterns, we have found the simulator to be quite useful in testing computer data processing techniques. The signal-to-noise ratio and the attendant statistical reliability of results obtained from complex modulated pulse trains can be determined by simulation for empirical estimates of the reliability of results obtained under similar conditions in the biological setting. We hope other investigators of neurological systems may find uses for this simulation technique.

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MAY 1976

266

ACKNOWLEDGMENT The authors are grateful to R. Disney for help with the statistical concepts and T. Birdsall for advice on the construction of the pseudorandom sequence generator. REFERENCES

[1] D. H. Perkel, G. L. Gerstein, and G. P. Moore, "Neuronal spike [2]

[3] [4]

[5] [6] [7]

[8]

[9] [101 [11]

[12]

trains and stochastic point processes-I: The single spike train," Biophys. J., vol. 7, pp. 391-418,1967. J. Tkvaril, T. Radil-Weiss, Z. Bohdaneck', and J. Syka, "Spontaneous discharge pattern of mesencephalic neurons: Interval histogram and mean interval relationship," Kybernetik, vol. 9, pp. 11-15,1971. C. K. Knox, "Signal transmission in random spike trains with applications to the statocyst neurons of the lobster," Kybernetik, vol. 7 , pp. 167-174, 1970. E. J. Bayly, "Spectral analysis of pulse frequency modulation in the nervous systems," IEEE Trans. Biomed. Eng., vol. BME-15, pp. 257-265, Oct. 1968. G. A. Korn, Random Process Simulation and Measurements. New York: McGraw-Hill, 1966. B. Elspas, "The theory of autonomous linear sequential networks," IRE Thans. Circuit Theory, vol. CT-6, pp. 45-60, Mar. 1959. A. J. Thomasian, The Structure of Probability Theory with Applications. New York: McGraw-Hill, 1969. A. Papoulis, Probability, Random Variables and Stochastic Processes. New York: McGraw-Hill, 1965. M. S. Bartlett, "Periodogram analysis and continuous spectra," Biometrika, vol. 37, pp. 1-16, 1950. - , "The spectral analysis of point processes," J. Royal Stat. Soc., vol. B25, pp. 264-296, 1963. W. J. Williams, S. L. BeMent, T. C. T. Yin, and W. D. McCall, Jr., "Nucleus gracilis responses to joint motion: A frequency response study," Brain Res., vol. 64, pp. 123-140, 1973. P. A. W. Lewis, Ed., Stochastic Point Processes: StatisticalAnalysis, Theory, and Applications. New York: Wiley-Interscience, 1972.

The Urinary Drop Spectrometer-An Electrooptical Instrument for Urological Analysis Based on the External Urine Stream ROGERS C. RITTER, NORMAN R. ZINNER, AND ARTHUR M. STERLING Abstract-This communication presents the design and an outline of clinical experience with two models of the urinary drop spectrometer (UDS), a noncontact optical instrument used for assessing urological abnormalities based on observation of the external urinary stream.

In previous publications the authors have discussed the early clinical applications and gross features of the urinary drop spectrometer (UDS) [ 1] -[5] . In general, a thin light sheet is sensed by a photodetector, and an electrical pulse is provided each time a drop in the external urinary stream passes through the light sheet. The drops have various sizes, shapes, rotations, vibrations, velocities and spacings, and the statistics of the optical observations of these properties constitute the analyzed data. These data represent the information extracted from the external urinary stream and have been shown, empirically, to be correlated with normal or various obstructed conditions of Manuscript received July 15, 1975; revised October 17, 1975. R. C. Ritter is with the Department of Physics, University of Virginia, Charlottesville, VA 22901. N. R. Zinner is with the Department of Urology, Charles R. Drew Medical School, Los Angeles, CA 90059. A. M. Sterling is with the Department of Urology, Academic Hospital, University of Leiden, Leiden, The Netherlands.

the urinary outlet tract. The accuracy and reliability with which the electronic signal represents the physical stream depends on design and construction features which have evolved over several years. In this communication the authors discuss the important aspects of this sequence of work, including a first, viable instrument, and a second, much improved spectrometer. I. BACKGROUND-MEDICAL INFORMATION Since urination is a highly volitional act [6] , one of the most important aspects of the UDS is that the instrument can be unobtrusive as well as physically noninvasive. Freedom from psychological interference is essential, and one goal is that the final instrument fit into a urinal and perhaps use infrared light. The normal adult male [5], [7] urinates from 100 to 800 cc at a voiding and has a mean flow rate from 10 to 40 cc/s. The mean volume of 250 cc is urinated on the average in about 12 s. Excluding small drops ("satellites" with volume < 0.02 cc) there are about 2100 drops having a mean volume of 0.12 cc (diameter = 0.61 cm). The average velocity at the meatus is about 200 cm/s, and the stream breaks up in a distance of 15 to 25 cm from the meatus. Changes in the drop ensemble, due to vibrations, collisions, etc., continue on downstream. In our experiments, these changes are essentially complete, as far as can be distinguished by time interval analysis of the drop train, at 60 cm from the meatus.

II. THE FIRST FIBER-OPTICS TRANSDUCERS The original transducers were designed around fiber-optics circle-to-slit converters in a configuration approximately as shown in Fig. 1. Several versions were used, one of them having a light field 12-cm wide and 25-cm long. Fibers were 0.01 0-in or 0.020-in diameter. Design of this system suffered most from complexity and unreliability. In spite of our best efforts, no manufacturer was found who could make the slits with uniformity better than about ±20 percent as measured by the output pulses from 1/4-in ball bearings (0.635-cm diameter) which, when dropped through the light fields, served as the standard occluding drops. Since fibers have a lobe-like intensity pattern, nonuniformity in fiber directions and spacings leads to light field nonuniformity. The parallel lobe patterns also lead to a fall-off at the edge of the light field, and mirrors at the sides provided an extension of the uniformity by the principle of imaging. Unfortunately, these front-surfaced mirrors were difficult to align and suffered also from mechanical and chemical instability. The photomultipliers, RCA 7102, were necessarily operated at several times their quiescent current ratings, which may have contributed to drift and noise. Drift, due also to decay in the intensity of the bulbs then used, was such that calibration checks and gain adjustments were needed twice a day. The use of photomultipliers in such an application is questionable, but simple phototubes gave such small pulses that we replaced them with photomultipliers. The nature of this problem is seen from the fact that a 0.12-cm3 sphere, of diameter 0.61 cm, occludes only 5 percent of the 12-cm width of the light field. In practice, the signal-to-noise ratio was at best 15/1 for such drops and it was more typically 12/1. The PM output pulses were a few tenths of a volt superposed on a 10-15 V dc level. Linearity was good, but ac coupling was needed to eliminate the dc component. This, in turn, made a base line clamp necessary. A factor present in any design of this instrument is the variation in ambient room lighting which finds its way into the receiving slits and causes changes in the output. In these early models there was insufficient collimation of received light, and room light from either tungsten or fluorescent lamps led to large 60 or 120 Hz baseline variations.

A flexible renewal process simulator for neural spike trains.

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, MAY 1976 262 1.0 ,....,... rz .? .0 (r&l -- ...- O.? r 700 .4.. J* -a 0-?1 0 .DI -do -03...
1MB Sizes 0 Downloads 0 Views