This preview shows pages 1–17. Sign up to view the full content.
This preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: I. 62 The Mysterious Change in the Time Constant: A Hypothetical Model The following is a very speculative systems block diagram proposal for how the brain changes
the time constant of the vestibular input. It's a little confusing, but it does involve known
anatomical wiring (i.e. for each location and pathway, there are know brain nuclei and fiber
tracts which could carry the signals which are hypothesized here), and does make use of an
important neural control principle in which a copy of the motor command signal is used
internally to anticipate subsequent changes in sensory inputs due to the subsequent movement. The following is the block diagram hypothesis. You are familiar with the head velocity input
(bottom) which is input to the semicircular canals (scc) with its transfer function with the canal
time constant (Tc=4 seconds for cats). This input goes to the vestibular nucleus (vn), is inverted O O
in sign and sent to muscles of the eye (the "plant"). The summation (E + H ) is the geometric
O
summation of head and eye velocity (or position) to get the gaze (G ) velocity (or position). Optokinetic System Proposed by Robinson not = nucleus of optic tract scc = semicircular canals
nrtp = nucleus regularis tegmenti pontis plant : ocular muscles and eye
vn = vestibular nucleus H = head position G = gaze position
W = world position of an object D = darkness L = light
e = retinal slip O
The novel part of the proposed system is that a copy of the signal E ’, the desired eye movement
velocity, is fed back as "velocity copy" to the nucleus regularis tegmenti pontis (nrtp), which
combines it with a signal coming from the retina which indicates how much the world (the visual
field) has slipped from its prior registration. Since the circuit involving the retina is slow, it is
useful to use velocity copy to anticipate the retinal slip due to the eye movement which is being carried out. l. 63 O . .
The net anticipated retinal slip (W h) is then low pass filtered and included as a dynamic Signal to
correct the oculomotor command. This part of the circuit is summarized below. nrtp eye vs. head
DOI + motion
3: ——————>
retinal plant
slip estimate of
~ world vs. world vs. retina motion head motion I
I
I low pass
I
I
I Hc T head velocity
estimate from semicircular canals We already know one of the parameters for this model: The time constant Tokan is 15 seconds.
(Remember the rotating drum experiment above, which is done in the light, but with the head still — H :0). We have one more measurement to make inside the drum, which will give us a measure of Gok. l.
2. Keep head still ———> no semicircular canal input. Keep frequency slow compared to oculomotor dynamics (Te)(w<6§ZIS—6C?; f<0.7 Hz). This makes the gain of the plant equal to 1.
1 15 secs Keep frequency slow compared to opto kinetic dynamics (Tokan) (0)< ; f<0.01 Hz) This means we can ignore the sTOkan term in the response. Result: E/W = 0.7 to 0.8; i.e., eye moves about 75% as much as the world moves. l. 64 Our simplified block diagram and calculation for GOk are as follows. 2.3— = 0.7 to O 8 W E _ Gok ' — l+GOk
anatomical W
evidence for Gok = 3 to 5
pathway Now we can calculate To, the time constant of the low pass filter, by manipulation of the transfer
functions in this block diagram. plant . l
Hok — l  k) ‘Gok
_ __ _ __ . . 9
E — [ To ] — STOkan + 1 (Can you derive this.)
s 1 _ k + l
but we know GOk 2 4 (3 to 5) therefore k z 0.75
also we know Tokan 2: 16 sec therefore T0 = 4 sec Finally, the mysterious change in time constant! One of the clever things this model does is to explain why the time constant of the VOR
changes from 4 seconds in the semicircular canals to 15 seconds in the vestibular nucleus and in O
the eye movement behavior. Let's calculate the transfer function of eye velocity (E' ) to head velocity as detected in the canals (lol C). £3.___l_______ sTO+1 16 4s+l TVOR STc+1
1‘31 lk TO ‘ 4 165+1 " TC sTVOR+l l. 65 MUSCLE SPINDLE MODELING Introduction The muscle spindles, along with the Golgi tendon organs, are sensors which provide information
to the central nervous system regarding the state of a muscle: its position (how far it has been
stretched), its velocity (how fast it is being stretched), and the tension being applied to it and
hence to the connections between it and the tendon. It follows that these sensors can provide
control information so that the animal or person has a good kinesthetic sense (knowing where a
limb is without visual information), good control over precise movements, and protection in case
of overload. From a control systems point of view it is useful to know how to characterize the
sensors in order to understand better their role in the feedback muscle control system. From the point of view of the course, we will use basic experimental data as an exercise in
system identiﬁcation, the process by which we choose the mathematical model and determine
the parameters which give it precision. From a linear systems modeling point of view, we intend that you learn three properties of these
as transducers: (a) the equivalence of "high pass filtering”, "differentiation", "phase advance" ~ all terms
used in linear systems theory — and "dynamic reSponding” or "phasic" (as opposed
to "static" or "tonic")  which are terms commonly used by biologists to describe
sensory responses. (b) how constrained the range of linearity may be in a real biological system (c) an appreciation for the muscle spindle's means of extending its dynamic range. You should refer to physiology texts for more information, including pictures like those shown
in lecture. I. 66 Physiological Overview
Figure reference: Kandel & Schwartz, 2nd edition, pp 451—456 (chapter 34) 1. Skeletal muscles are connected by tendons to bone. The Golgi tendon organs sense tension
and can be considered force sensors. 2. Muscles contain extrafusal fibers, which contract and cause motion of the limb. These are
innervated by "alpha" motor neurons, which have the largest and fastest conducting axons in the nerve bundle which innnervates the muscle. 3. The intrafusal fibers, mixed in with the extrafusal fibers, contain the muscle spindles which
respond to motion or change in position. In addition these fibers can contract in order to change the sensitivity of the spindles. 4a. The nuclear chain fibers are "static" responding: ideally, they fire proportionally to the
position or length of the stretch. 4b. The nuclear bag fibers are "dynamic" responding: ideally, they fire proportionally to the
velocity at which the muscle is being stretched. 5a. The static gamma motor fibers innervate the nuclear chain fibers, causing them to contract,
changing their senstivitiy. Example: when the motor fibers fire at the correct rate, the sensor
reports small changes in position with small changes in firing rate; when the motor fibers are
shut down experimentally, the intrafusal fiber goes slack and the sensory fibers do not fire at all,
as they don't change length, even though the extrafusal muscle does change length. 5b. The dynamic gamma motor fibers innervate the nuclear bag fibers, causing them to contract,
changing their senstivitiy. 6. Unfortunately, there are limits to our ability to make recordings in which we can distinguish
among the different motor and sensory fibers. There are two principles which can help, however: (a) anatomically, it is sometimes possible to separate groups of nerve fibers within a
bundle. Unfortunately, most large skeletal muscles (e.g. gastrocnemius) are mixed
nerves where this is somewhere between difficult and impossible. (b) electrically, large axons have both bigger signals and faster conduction speeds. It is
difficult to take advantage of the size of the signals, because one does not know ahead of
time how many neurons in each class are responding (recorded size=(# neurons)*(size of
signal from each). It is often possible to distinguish by conduction velocity, by placing
two electrodes along the nerve. Although it is possible to distinguish three partially overlapping classes of afferent (e.g. sensory)
muscle fiber signals (primary spindle afferents, secondary spindle afferents, and Golgi tndon
organ afferents), it is not possible to separately record signals from dynamic versus static
responding fibers. It is quite likely that the nervous systems knows which ones are dynamic and
static, but the control systems engineer, who might wish to know the activity of these neurons in
order to control a prosthetic device, cannot find them directly. I. 67 Modeling/ System Identiﬁcation Question: a priori, can you suggest a simple transfer function model which relates neural firing
rate to changes in input postion? I.e., what is the Laplace transform of the operator which
defines the ideal properties of the muscle spindle transducer? An equivalent question is, do you
expect the spindles to act as low pass filters, high pass filters, ..,? A reasonable, but wrong answer, is that they act as low pass filters, since no response can occur
instantly. This may be true, but only in that it is possible for an experimentalist to drive the
muscles at frequencies which are beyond the physiological range. The key to the answer is that the response is proportional to the sum of a set of fibers which are
static and a set which respond dynamically. rate = {static responders} + {dynamic responders}
r(t) = Ka s(t) + Kb d(t) The input is an experimentally controlled function x(t) which we hope mimics real behavior
where the muscle is controlled by the animal’s own alpha motor fibers. Hopefully: K1 s(t) = Ks x(t) and K2 d(t) = %
and r(t) = r(t) 2 KS x(t) + 19%
so that
_ E“; _ _S_
G(s) — x(s) — K (1 + CDC) Ref: Matthews&Stein. Sensitivity of muscle spindle afferents. J.Physiol.(l969),ﬂ,723745. These figures show temporal data and Bode plots relating muscle spindle firing rates to a
stimulus which changes the length of a muscle. Figure 1 illustrates how the data appear: average firing rate versus time for a sinusoidal stimulus
with frequency of 1 Hz and amplitude of l millimeter. Neural firing rates (impulses / second)
vary from 10 to 60, and roughly follow the stimulus. They show phase advance —— i.e. the peak
in the firing rate occurs roughly 45 degrees before the peak muscle extension. This is a natural
consequence of the physiology: static responders have a peak response at maximal muscle
extension; the dynamic responders have a maximal response at peak muscle velocity, which
occurs 90 degrees earlier. Since the primary afferent nerve includes both, the average phase
response is in between, and depends on the relative amplitudes of the groups of fibers. The next figure shows a Bode amplitude plot for both the primary and secondary spindle
afferents. Each shows a response which is constant in amplitude up to approximately 1.5 Hz,
and linearly increasing thereafter. Further, the slope of the plot is 20 dB/decade, or a factor of 10
increase in firing rate for each factor of ten change in frequency. Each can be modeled by the transfer function above, with 0)C=21tfc=10 rad/sec, and values of K equal to 100 impulses/sec per
millimeter for the primary afferents, and 8 impulses/sec per mm for the secondary afferents. R
G(s)=X%=K(1+$sC) l. 68 System Identiﬁcation:
We have already done system identification for the muscle spindles. The basic steps were: (a) determine that the system can be modeled as a linear system over some range of input and
output behavior (the authors have done this) (b) determine the form of the system transfer function:
(i) by first principles, as argued from the observation that the nerve consists of a mix of dynamic and static responders
(ii) by inspection of the Bode plot, in which we identify the slope of the gain curve over different frequency ranges (c) use the Bode plot to calculate the parameters K and (DC. Limitations to the Applicability of the Model
The article provides several views of the limitations of this linear model: (a) casual observation of the spike rate vs. time function shows that there is considerable
variation from linearity. (b) in the Bode amplitude plot, there is some deviation from the model, especially at high
frequencies. (However, think of what it means to move a muscle at 100 to 300 Hz!) (c) in the Bode phase plot, the deviations are much worse (d) the response amplitude is linear with stretch amplitude up to 2 mm for secondary fibers and
up to 0.15 mm for primary fibers, after which the gain falls off. Sensitivity Changes with Gamma Stimulation Figures taken from: Chen and Poppele — Small signal analysis muscle spindles with fusimotor
stimulation, J. Neurophysiol. (1978), 5L1, 15—27. Several graphs are shown which document how the transfer functions change when either the
dynamic or the static efferent (i.e., motor) fibers are stimulated. When the static gamma motor
fibers are stimulated, the sensitivity decreases; when the dynamic gamma motor neurons are
stimulated, the sensitivity decreases. These graphs argue that, once the gamma motor input is
fixed, the system is linear. However, the dependence on the gamma motor input is hard to model as linear. Further, the gain is a complicated function of the peak to peak stretch amplitude. I. 69 LINEAR, N ON RATION AL SYSTEM This example of a linear model of a biological system is that of a non—rational system. It
is so named because the exponent of the term in the Laplace transform is not an integer,
resulting in a nonrational operation on the variable. For instance the value may be l/2, which indicates that the square root is being taken. Such systems arise naturally, perhaps much more than we acknowledge, as the non—rational form
indicates generally that a system has a distribution of values, rather than a single value. For
instance, assume that you have a collection of low pass filters operating in parallel such that the
output is the sum of all the filtered inputs. Assume further that the DC gains are the same, but
that the cutoff frequencies vary from 1 Hz to 10,000 Hz, with an equal distribution over this
range. Assume you create a Bode plot of the output vs. the input. You would find that the
output is constant from 0 to 1 Hz, and that it slopes at —20 dB/decade beyond 100 Hz, with a
phase shift of —90 degrees. In between 1 Hz and 10,000 Hz, you would find that it would slope at
—10 dB/decade, with a phase shift of —45 degrees. Over the range from 0 to 10,000 Hz you might
characterize it as having a transfer function of: l \Il+s The biological example here is a model of the slowly adapting crayfish stretch receptor, and is
taken from the work of MC. Brown and R. B. Stein, J. Cybernetics, 1966. This same approach
has recently been used by Prof. Anastasio of the Physiology Department at the UIUC for modeling parts of the VOR. or (l+s)1/2 The basic experimental setup is illustrated here, where a crayfish muscle is isolated from the
crayfish and stretched under careful control, while muscle tension, length, and neural firing rate were measured. S train
Gauge Muscle
Tension Muscle Optical
Puller Length
Measurement NerVC /'
l Output
Input Amp Rate I. 70 Experiment 1 The responses to a step change in length are illustrated here. The first goal was to
determine whether the nerve fibers were responding to changes in length or to tension. The
researchers found a strong correlation to both. making it impossible to distinguish between the two hypotheses. On the other hand, both responses were highly linear. This may seem odd at first, since the pulse
rates decline in time, but it was true that double the step change in length led to twice the firing
rate over the entire time course of the response. System Appears Linear X2 stretch
0.5 m1 X1
tension
correlation
94%
correlation
95%
pulses
per
second . _
What is receptor transduc1ng?
land Length or tension? Replot Data The experimenters then replotted the data on a log/log scale. This is a classic approach to
discovering power law relationships, which you hopefully learned in physics or in mathematics. linear on loglog scale
100 / => power law relationship k l— —l— time (log scale) .1 1.0
sec sec firing rate
log scale
N L11
0 o .._.
O I. 71 The firing rate curve now appeared linear on the log/log scale, indicating that the firing rate
equation could be written as: y(t) = a & t'k or, (pulses per second) = (constant) (length of step input) (time)(fmcti0“) The function y(t) actually has a Laplace transform” Y(s) = a Q L[t’k] (where L[ ] is the Laplace operator) Y(s) = 3 9r l"(l—k) 51"] (where l" is the Gamma function; it is available
in Matlab; for example, T(.79)=1.l757) Now, the input
x(t) =
X(S) = A
u(t) x = the length of the step, a constant for any trial /s ><>><> And the system transfer function G(s) is given by A _ k~1 =X(_S):£Z‘_£L_k)s_ : ar(1_k)sk = Csk A
El constant
G(s) = C sk , where k is not necessarily an integer The following are experimental values of the terms in the transfer function as determined from
step response Values of k Values of a
(a dimensionless constant) (impulses/sec per mm of stretch)
range 0.157 —> .244 20.2 —> 62.6
mean 0.210 41.9 std. error i0.015 i7.3 I. 72 Frequency Response and Ramp Response: As stated earlier, for a linear system, if you know the step input, then you can compute the
impulse response and the frequency response, as well as the response to any other input. What
does the frequency response of this nonrational system look like? Frequency Response: G00» GOOD“ = C wk (Dk
= c wk(ej90°)k = C madam For some reason, many students have a hard time recognizing that j = \fT = elW/Z = l exp(90°) (a unit vector pointing "up")
(j)k = 1 exp( k 900 ) (a unit vector point to wherever k times 90° is) Thus, for the stretch receptor data, where k = 0.21
lGQw)l = C (0021 q) = 905021 = 18.9° i.e., for x(t)= Acos wt
y(t) 2 [AC 03021] cos [(0t+18.9°] Prediction phase
constant 189° Log
Gain log f This is in agreement with our previous model of the muscle spindle — this is a high pass filter
consisting of a mixture of static and dynamic responding fibers. However, here the cutoff
frequencies of the individual fibers varies over a wide range  which is actually a much more
realistic situation than in the previous muscle spindle model. l. 73 Ramp Input Prediction V The transform of a ramp is —2 where V is the velocity of the ramp. Thus the response should
s be: G(s)  V/52 The experimental input was a ramp of duration t1, which can be modeled as the difference
between two ramp functions, one of which has a pure time delay from the other, as indicated in
the drawing. ramp

K slope velocity = —
= V = x/tl
t1 Thus, the input transform is: X(s) = {X2 — X2 (e’5t1))
s s X(t) 71> VCsk
52 Va F(1—k) _[ V ,.,_
=—';EI—(1—e“)= cw Va Va
_ 11 _ lk
y“) 11 t “ l—kt Y(S) = X(S) G(S) = (l — e—S‘I) ; (remember, C = a F (l—k) ) use F(n +1): n Rn) tStl y(t) 1_k [tlk—(t—t1)1k] t>t1 Prediction 1: plot log (y(t)) vs. log t, slope = l — k i . Vatll'k aQI—k
Prediction2: maxy(t)=y(t1)=—1—:—k—' = 1_k vk plot log (max(y(t))) vs. log (V) I. 74 Crayfish Photoreceptor Transfer Function: a System Identification Example
Reference: Stark, Neurological Control Systems, Chapter 1. Introduction: The lowly crayfish has a photosensitive ganglion located near its tail (this was shown by Ladd
Prosser, former head of the Physiology Department at the UIUC). They are used to help the
animal avoid light and to take cover so as to avoid predators. (Limulus, the horseshoe crab, also
has photosensors along its tail region, which seem to have sufficent resolution to allow them to
recognize the presence or absence of limulus sized objects within a few feet, in tide waters under full moons when they mate!) One approach to modeling their response has been to focus a light on the photosensitive
ganglion, modulate the light with sinusoids, step and ramp functions, and to measure the spike
rates of the abdominal nerve which connects the ganglion to the animal's head. The pulse rates
are shown to be very roughly linear, but with a twist not usually encountered in undergraduate
systems courses: there is a pure time delay, which can be represented exactly using Laplace transforms. Our use of this example is for system identification with a time delay. It is an example of a Non
Minimum Phase (NMP) system. In a Minimum Phase system, once you know the amplitude
Bode plot, you can write down the transfer function, and you know what the phase plot looks
like. In an NMP system, the gain (amplitude) is the same, but there is an additional phase shift.
The only source of such a shift considered in this course is a time delay. There are many others,
and they are often used in electronic filter design to obtain desired filter properties. Instrumentation Although the figures highlighting the signals attained during these experiments are old (30
years!), the essentials of the processing are not. Now, of course, everything would be done in
software, whereas the neurophysiologists of old had to be more clever. The instrumentation is
able to distinguish the response of (a) a population of neurons  all neurons whose spikes have an amplitude above a
minimum (b) a single unit —— a single (hopefully) neuron whose spike amplitude falls in a narrow
range I. 75 Step Response A linear system can be characterized by its step response. Graphs are shown which indicate that
a step change to 5 millilumens/mm2 is too great for linearity to hold but, a step change to 2 millilumens/mm2 is within linear region. There appears to be a delay of between 1 and 2
seconds in the onset of the response to the step input. Be mindful of the following Laplace relationships: A step input, u(t) has the transform Si . Let the step response be h(t), with transform H(s).
Let the impulse response be g(t), with transform G(s).
Then: h(t) = g(t) * u(t) (convolution)
H(s) = G(s) Si
and
G(s) = s H(s) Thus, to use this method to find the system transfer function, we have to characterize the time
domain signal h(t), take its transform, multiply by s to get G(s), and, if we want the impulse
response, inverse transform to get g(t). I. 76 Sinusoidal Response Population Response: Graphs are shown indicating the response to sinusoids. Bode amplitude and phase plots are
shown for the population response which indicate a first order, low pass filter function with a
break frequency of 0.14 Hz (0.9 radians/sec). Below are Matlab generated versions of the Bode
plots, with and without the extra phase shift. Note that, although the plots here include the gain
363 = 50 dB, often plots are normalized so that the maximum gain is 1 or 0 dB. 10" 10° 10‘
Frequency (HZ) —3BD  L (A  i 17';
10‘3 1.0" 10" 10° 10‘
Frequency (Hz) Bode plot with both: Minimum Phase only
(upper phase plot), and with NonMinimum
Phase term (lower phase plot). . 363
G(s):1+1.1s 363 e'O55
i+i.iS Minimum Phase: Total Transfer function: G(s) = The transfer function is a ratio of pulse rate to
light intensity: P(s) _ ‘ _ 363 e'0~55
m ‘G“)‘ (1+1.1 s) 55r—— ALA
mo Galn In (18 1:.
o i T.  i x l l l
O 2 0 4 0.6 08 1 1 2 14 1 6
Frequency In Hz Phase in Degrees _J..__. .__._l __l ‘ l l
0 B 0.8 1 1,2 14 1 6
Frequency m Hz Total transfer function plotted against a linear
frequency axis.
—0.55
(3(8) = 363 e l+l.ls Refer to earlier notes on how to calcuate the
time delay from this Bode plot. The time domain function relates pulses per
second, P(t), to the light intensity L(t) in millilumens per mmz: 363 L(t—0.5) = 1.1 Cu) 3+1) 300
200p ﬂ,»‘”" i”“*as_ J 1007/ lmag AXIS
o 100— 200   l
100 200 300 400
Realexs op. ‘aqgoo Single Unit Response: V 10‘ 10c 10” 10°
Frequency (Hz) to“ 10Q 10* 10° 10
Frequency (Hz) 1. 77 Often, the Bode plot data are replotted as a
Nyquist plot, where the distance from the
origin is the amplitude, the angle is the phase
angle, and the parameter is frequency. For
instance, the value at 0 frequency is 363 with a
0 degree phase angle. At 0.2 Hz, the phase
angle is 90 degrees and the amplitude is 220 or 47 dB. The first part of the curve (mostly in the
bottom half) represents evaluating the transfer function for positive frequencies jot), while the
second part (mostly the top half) represents evaluation of negative frequencies j(D. A similar analysis was done for the single unit
data, although the plot indicated that the slope
was 40 dB per decade and that the system was
a second order system. There is a time delay
of 1.02 seconds. The transfer and time domain functions are: _ _ 3261.05
Us) "G(S)_ (1+1.2 s)2 d2P dP
32 L(tl)— 1445;2— +2.4 E +P I. 78 Using both Sinusoidal and Step Response Data If you assume that a system is linear, then:
if you know the frequency response, you can calculate the step reponse (and the impulse
response)
if you know the step response, you can calculate the frequency respsonse So why do both?
It is a way to check to see if the system really is linear. If you check the responses in
Stark's article, you'll be skeptical about the time delay component. Transfer Function Impulse Response Transform of the Step Step Response
G(s) Response (H(s) (h(t))
(The transform of the
impulse response) #1 +1 I1.1s em s<1+11.1s> l'em']
1.1 e'O55 e—(t—0.5)/1.l 3:955— 1_e—(t—0.5)/l.l
l+1.1s S(l+115)
6.105 t e(tl.0)/1.2 el ~05 1_e—(t—1.0)/l .2
(1+ 1.2s)2 s(1+1.2s)2 _ ﬁemoyrz Linearizing a System about its Operating Point Modelers of biological systems use a traditional technique, which is to linearize a system about
its operating or resting point. The general argument is that if y is a function of x, then a Taylor
series expansion can be made about any given operating point. y {linear region about x0
y / Taylor expansion
0 '(X—XO)+...
x x0
x0 hopefully
ignorable This has arisen in the models of the neuronal firing rates of the VCR, where it is usually assumed
that every neuron has a sponstaneous, steady state firing rate which is modulated up or down. In the crayfish model a step change of 5 millilumens/mm2 is too great for linearity to hold, but,a
step change of 2 millilumens/m2 is within linear region. In the modeling done here, the chosen
operating point where the light level is 2.5 millilumen/mm2 and the resulting operating range is
[0.5 millilumen/mmz, 4.5 millilumen/mm2 ). ...
View
Full
Document
This note was uploaded on 01/23/2012 for the course EEL 6502 taught by Professor Principe during the Spring '08 term at University of Florida.
 Spring '08
 PRINCIPE
 Signal Processing, The Land

Click to edit the document details