EP Model -- Delphi version, revised

[From Bruce Abbott (2014.02.05.1130 EST)]

I’ve resolved the issue about the wrong muscle being contracted to resist a gravitational load. It looks like the original paper had reversed the order of angle subtractions in two equations and then reversed the subtraction in the final torque computation, which corrected the direction in which the modeled joint moves. It’s similar to computing the error, e, in a control system by subtracting r from p instead of the reverse and then reversing the sign of the output to compensate.

I also realized that the Lan and Zhu (2007) paper was not based on Asatryan and Feldman (1965) as I had misremembered but rather on St.-Onge, Qi, and Feldman (1993); I’ve corrected that in the main form’s label. The corrected version of my program has been substituted for the original on my Google website (https://sites.google.com/site/perceptualcontroldemos/home/other-demos ), in case you want to download the new version.

Bruce

[From Rick Marken (2014.02.05.1515)]

Bruce Abbott (2014.02.05.1130 EST)–

BA: I’ve resolved the issue about the wrong muscle being contracted to resist a
gravitational load.

I think the prior version was a better replication of the Lan/Zho EP model. The main reason is that the new model doesn’t seem to work correctly when

C is greater than about 10. And neither version of the model seems to work as well as your spreadsheet version. The spreadsheet version seemed to behave just like the Lan/Zho model for the combinations of R and C values that they reported and that you tested. I show this below: Hope the graphic comes through:

The top graph is the Lan/Zho plot of time variations in limb angle for various settings of R and C. The graph below is your graph of the same thing from your spreadsheet model for the same values of R with C set to 5. The match seems perfect. The results for you spreadsheet model with R = 60 and when C= 80 (which gave the most direct approach to the new angle with no wobble) was also a match to the Lan/Zho data. Your first Delphi version also seemed to behave this way; it gave the best results (in terms of stability) with C=80 for all values of R. Unlike the Lan/Zho model, the new Delphi version of the model doesn’t work at all with C=80 for all values of R.

One problem with both versions of the Delphi model is that the limb angle axis runs from -90 to 90 while in the Lan/Zho and spreadsheet graphs angle runs from 0 to 140. Maybe that’s part of the problem with the new model.

But I think you were on the right track with the prior version of the Delphi model in the sense that it seemed to behave exactly like the Lan/Zho model with their settings of R and C. I think the correspondence between the behavior of your model and theirs should be the basis for determining whether or not your model is equivalent to theirs.

I agree that you really have to get this right before you can say anything about whether the EP model controls or not. Right now neither version of the Delphi model controls so I’m pretty sure that the EP model is not a control model. But either way, it’a a win for PCT; if the EP model doesn’t control we can show easily that it can’t possibly be a model of limb movement and if it does control we can point out that it does so because it perceives limb position; so it’s a PCT model after all, just with different initials;-)

So keep up the good work.

Best

Rick

···

Richard S. Marken PhD
www.mindreadings.com

The only thing that will redeem mankind is cooperation.
– Bertrand Russell

[Rick Marken (2014.02.05.1520)]

The graph didn’t show up in my copy so I’m attaching it here, just in case.

Best

Rick

···

On Wed, Feb 5, 2014 at 3:15 PM, Richard Marken rsmarken@gmail.com wrote:

[From Rick Marken (2014.02.05.1515)]

Bruce Abbott (2014.02.05.1130 EST)–

BA: I’ve resolved the issue about the wrong muscle being contracted to resist a
gravitational load.

I think the prior version was a better replication of the Lan/Zho EP model. The main reason is that the new model doesn’t seem to work correctly when C is greater than about 10. And neither version of the model seems to work as well as your spreadsheet version. The spreadsheet version seemed to behave just like the Lan/Zho model for the combinations of R and C values that they reported and that you tested. I show this below: Hope the graphic comes through:

The top graph is the Lan/Zho plot of time variations in limb angle for various settings of R and C. The graph below is your graph of the same thing from your spreadsheet model for the same values of R with C set to 5. The match seems perfect. The results for you spreadsheet model with R = 60 and when C= 80 (which gave the most direct approach to the new angle with no wobble) was also a match to the Lan/Zho data. Your first Delphi version also seemed to behave this way; it gave the best results (in terms of stability) with C=80 for all values of R. Unlike the Lan/Zho model, the new Delphi version of the model doesn’t work at all with C=80 for all values of R.

One problem with both versions of the Delphi model is that the limb angle axis runs from -90 to 90 while in the Lan/Zho and spreadsheet graphs angle runs from 0 to 140. Maybe that’s part of the problem with the new model.

But I think you were on the right track with the prior version of the Delphi model in the sense that it seemed to behave exactly like the Lan/Zho model with their settings of R and C. I think the correspondence between the behavior of your model and theirs should be the basis for determining whether or not your model is equivalent to theirs.

I agree that you really have to get this right before you can say anything about whether the EP model controls or not. Right now neither version of the Delphi model controls so I’m pretty sure that the EP model is not a control model. But either way, it’a a win for PCT; if the EP model doesn’t control we can show easily that it can’t possibly be a model of limb movement and if it does control we can point out that it does so because it perceives limb position; so it’s a PCT model after all, just with different initials;-)

So keep up the good work.

Best

Rick


Richard S. Marken PhD
www.mindreadings.com

The only thing that will redeem mankind is cooperation.
– Bertrand Russell


Richard S. Marken PhD
www.mindreadings.com
The only thing that will redeem mankind is cooperation.

                                               -- Bertrand Russell

[From Bruce Abbott (2014.02.05.1835 EST)]

Rick Marken (2014.02.05.1515) –

Bruce Abbott (2014.02.05.1130 EST)

BA: I’ve resolved the issue about the wrong muscle being contracted to resist a
gravitational load.

RM: I think the prior version was a better replication of the Lan/Zho EP model. The main reason is that the new model doesn’t seem to work correctly

BA: You may be working with the wrong version. After posting the revision, I discovered that it wasn’t handling the added weight correctly. I made another change which seems to have fixed that problem.

Bruce

[From Rick Marken (2014.02.05.1830PST)]

Bruce Abbott (2014.02.05.1835 EST)–

RM: I think the prior version was a better replication of the Lan/Zho EP model. The main reason is that the new model doesn’t seem to work correctly

BA: You may be working with the wrong version. After posting the revision, I discovered that it wasn’t handling the added weight correctly. I made another change which seems to have fixed that problem.

RM: It sure did! It works great and it does indeed seem to control joint angle. I would like to know what variable is controlled by the model. Is it the angle itself or some correlate of the angle? The R parameter seems to function like a reference for angle so there must be somewhere in your code where R is compared to that angle (or its correlate). And C seems to function as a gain parameter. The higher C, the tighter the control. We could probably figure out how the control loop works from the code but it might be easier if you could posy a functional flow diagram of the model. I’d really like to see how it works.

Best

Rick

[From Bruce Abbott (2013.02.06.1015 EST)]

Rick Marken (2014.02.05.1830PST) –

Bruce Abbott (2014.02.05.1835 EST)

RM: I think the prior version was a better replication of the Lan/Zho EP model. The main reason is that the new model doesn’t seem to work correctly

BA: You may be working with the wrong version. After posting the revision, I discovered that it wasn’t handling the added weight correctly. I made another change which seems to have fixed that problem.

RM: It sure did! It works great and it does indeed seem to control joint angle. I would like to know what variable is controlled by the model. Is it the angle itself or some correlate of the angle? The R parameter seems to function like a reference for angle so there must be somewhere in your code where R is compared to that angle (or its correlate). And C seems to function as a gain parameter. The higher C, the tighter the control. We could probably figure out how the control loop works from the code but it might be easier if you could posy a functional flow diagram of the model. I’d really like to see how it works.

BA: Like our own models, this version of the EP model leaves out a lot of detail that a complete model will need to include. For example, the muscles are not explicitly modeled; instead, functions represent the joint angle that would be realized if the muscle attained a particular length without constraint from the opposing muscle. So this version of the EP model uses these joint angles as proxies for muscle lengths. The reference muscle lengths for the flexor and extensor muscles can differ even though there can be only one actual joint angle at the joint operated by those muscles, reflecting different “commanded” lengths. The muscle is assumed to contract in an effort to move the joint to the specified angle, but counter-torques from the opposing muscle or outside disturbances may prevent the specified angle from being reached. The result is that at equilibrium the two muscles may be co-contracted and producing opposing but balanced torques on the joint; consequently the joint remains stationary.

BA: In the model there are two “commands” sent from elsewhere. “R” represents the commanded joint angle and “C” the degree of co-contraction of the opposing muscles. These two command variables are combined to produce commanded values of a variable called “lambda” associated with each muscle. For one muscle, lambda is the sum of R and C, for the other, lambda is their difference. Lambda represents the threshold values of muscle length at which alpha motor neurons that operate the muscle begin to produce impulses that contract the muscle. As this threshold is exceeded, more and more alpha motor neurons begin to fire impulses and the frequency of impulses also increases.

BA: The ACTUAL values of lambda (what the authors call its “dynamic” values) are functions of the “commanded” values but are also modified by the rate of muscle contraction and reciprocal inhibition (via inhibitory interneurons in the spinal cord) arising from the activity of the opposing muscle’s motor neurons. The influence of the rate of contraction is modeled as a damping factor that is proportional to the rate of change of joint angle. (In the model this is a constant of proportionality, mu, multiplied by joint velocity, omega.) Reciprocal inhibition is represented by the variable rho, which is computed in the model before being applied to the computation of the dynamic lambda. The computations of rho are similar to those of dynamic lambda for the opposing muscle but are based on the difference between the commanded lambda for that muscle (in degrees of joint angle) and the actual angle, theta. Also included are the same damping factor as before and a variable called “h.” H is in degrees and is supposed to represent a difference in the activation thresholds of interneurons, relative to motor neurons. (In the simulation it is a constant +10 degrees). If the computed rho is negative, it is set to zero in the model, because neurons cannot have a negative firing rate.

BA: The dynamic lambda depends on both the commanded lambda and the rate of change in joint angle (due to damping, which is proportional to rate), as well as reciprocal inhibition. This dynamic lambda represents the dynamic thresholds for motor neuron activation. Muscle activation is computed from the difference between the actual length of the muscle (represented by actual joint angle) and dynamic lambda, multiplied by a constant of proportionality. If the value is negative, it is set to zero (because there can be no negative muscle activation). If it is positive, it becomes the exponent of an exponential function. Thus muscle activation is modeled as an exponential function of the difference between dynamic lambda and the actual joint angle. This computation is carried out separately for the flexor and extensor muscles. The “actual joint angle” as well as some of the difference between the static and dynamic values of lambda presumable arises from the output of the muscle spindles, sensors in the muscle sensitive to muscle length and rate of change in length. The difference between the dynamic activation threshold (represented as a joint angle) and the actual joint angle is the error that is corrected in the model via negative feedback – muscles are activated to produce torques in a direction that reduces the difference.

BA: Muscle activation, E, is multiplied by a constant of proportionality to produce the torque generated on the joint, which acts in opposing directions for the flexor and extensor muscles. The difference in torques generates an angular acceleration of the joint. This acceleration depends on the net torque, divided by I (that’s “i” not “one”), the moment of inertia of the forearm. (This is the rotational analog to Newton’s a = f/m (i.e., acceleration equals force divided by mass).

BA: Joint rotation accelerates under this net torque, this generates a change in velocity at each time-step of the simulation dt, producing a new velocity, and the change in velocity produces a new position of the joint. The joint rotates until the torques balance. In the absence of an external disturbance or encountering the limits of joint motion, the forearm moves to the angular position specified by R and maintains a level of co-contraction determined by C.

BA: I hope my explanation is clear enough without presenting the actual formulae. These are available in the Lan and Zhu (2007) paper. The values of all constants in the simulation were also taken from that paper, page 200.

Bruce

[From Rick Marken (2014.02.06.1045)]

Bruce Abbott (2013.02.06.1015 EST)--

RM: We could probably figure out how the control loop works from the code but it
might be easier if you could posy a functional flow diagram of the model.
I'd really like to see how it works.

BA: Like our own models, this version of the EP model leaves out a lot of
detail that a complete model will need to include.

RM: Yes, but unlike our own models it is very difficult to understand
how it works. There were inconsistencies between the Lan/Zho equations
and diagram that were quite puzzling and exasperating. Your first go
at a model, presumably based on the Lan/Zho equations, produced
response curves just like those seen in the Lan/Zho paper; so it
looked like you nailed their model. And when we tested that model for
control (in both the spreadsheet and Delphi versions) by adding a
weight, there was no control.

Now you re-write the model based, apparently, on another paper (by St.
Onge et al) and suddenly it controls. While the response curve
produced by the new model seem similar to the ones shown in Lan/Zho,
the behavior of the variables in the model is quite different that in
the Lan/Zho version of the model. In particular, teh Lambda and
Lambda* values behave completely differently in the two versions of
the model. In the Lan/Zho version (which you call Asatryan/Feldman in
the Delphi version) Lambda* seems to be proportional to the
difference between rho and Lambda, suggesting the rho is the
controlled variable (which it looks like in the Lan/Zho diagram) and
Lambda* is an error. In the St. Onge et al version it's Lambda* that
seems to be the controlled variable and rho the error. So why that
model controls and the Asatryan/Feldman models doesn't is a puzzle.

>BA: For example, the muscles

are not explicitly modeled; instead, functions represent the joint angle
that would be realized if the muscle attained a particular length without
constraint from the opposing muscle. So this version of the EP model uses
these joint angles as proxies for muscle lengths. The reference muscle
lengths for the flexor and extensor muscles can differ even though there can
be only one actual joint angle at the joint operated by those muscles,
reflecting different "commanded" lengths. The muscle is assumed to contract
in an effort to move the joint to the specified angle, but counter-torques
from the opposing muscle or outside disturbances may prevent the specified
angle from being reached. The result is that at equilibrium the two muscles
may be co-contracted and producing opposing but balanced torques on the
joint; consequently the joint remains stationary.

RM: This doesn't help me understand how the control is achieved in the
St. Onge version of the model. And there definitely is control going
on; when you increase the downward torque on the limb (by increasing
weight) the flexor increases upward torque to nearly exactly
compensate for the downward torque; and it you increase the upward
torque (by making the weight negative; equivalent to pushing up on the
limb) the extensor increases it's downward force to compensate. There
is control going on and what is controlled is clearly Lambda*; that
variable is being brought into a match with Lambda. Rho now functions
as the error signal. This is unquestionably an active control process
since the normal response to a downward torque on a limb on a hinged
joint maintained at a particular angle is not increased but decreased
tension in the flexor spring. So I think it's still important to
figure out why this version of the St. Onge version of the model
controls and why the Lan/Zho version doesn't. Actually, I think what
we still need to figure out is what is the EP model? Actually, read
on; I think I figured it out!!

BA: In the model there are two "commands" sent from elsewhere. "R"
represents the commanded joint angle and "C" the degree of co-contraction of
the opposing muscles. These two command variables are combined to produce
commanded values of a variable called "lambda" associated with each muscle.
For one muscle, lambda is the sum of R and C, for the other, lambda is their
difference. Lambda represents the threshold values of muscle length at
which alpha motor neurons that operate the muscle begin to produce impulses
that contract the muscle. As this threshold is exceeded, more and more
alpha motor neurons begin to fire impulses and the frequency of impulses
also increases.

RM: Yes, lambda does seem to function like a reference signal; it
seems to "command" the value of lambda* (at least in this new version
of the model; it seems to command rho in the earlier version).

BA: The ACTUAL values of lambda (what the authors call its "dynamic" values)
are functions of the "commanded" values but are also modified by the rate of
muscle contraction and reciprocal inhibition (via inhibitory interneurons in
the spinal cord) arising from the activity of the opposing muscle's motor
neurons. The influence of the rate of contraction is modeled as a damping
factor that is proportional to the rate of change of joint angle. (In the
model this is a constant of proportionality, mu, multiplied by joint
velocity, omega.) Reciprocal inhibition is represented by the variable rho,
which is computed in the model before being applied to the computation of
the dynamic lambda. The computations of rho are similar to those of dynamic
lambda for the opposing muscle but are based on the difference between the
commanded lambda for that muscle (in degrees of joint angle) and the actual
angle, theta.

RM: Ah, there it is! Theta is the controlled variable. Somehow theta
must be "allocated" to lambdaF* and lambdaE*, the muscle length
perceptions. So theta is not controlled at all; what is controlled is
the relative muscle lengths, lambdaF* and lambdaE*. The references for
these muscle lengths, lambdaF and lambdaE, which are set by R and C,
are calibrated so that a particular combination of these reference
equal a particular limb angle. It looks like a particular ratio of
lambdaF to lambdaE corresponds to a particular angle, but it may be a
difference.

Anyway, that's what's going on in the St. Onge version of the model.
It's a control model controlling perceptions of flexor and extensor
muscle length and, as a side effect, controlling objective limb angle.
So EP is just an incoherent version of PCT with an attitude;-)

Again, I think this is a great basis for a paper!! Feldman, in railing
against PCT, was just railing against his own model. I love it!

Best

Rick

Also included are the same damping factor as before and a

···

variable called "h." H is in degrees and is supposed to represent a
difference in the activation thresholds of interneurons, relative to motor
neurons. (In the simulation it is a constant +10 degrees). If the computed
rho is negative, it is set to zero in the model, because neurons cannot have
a negative firing rate.

BA: The dynamic lambda depends on both the commanded lambda and the rate of
change in joint angle (due to damping, which is proportional to rate), as
well as reciprocal inhibition. This dynamic lambda represents the dynamic
thresholds for motor neuron activation. Muscle activation is computed from
the difference between the actual length of the muscle (represented by
actual joint angle) and dynamic lambda, multiplied by a constant of
proportionality. If the value is negative, it is set to zero (because there
can be no negative muscle activation). If it is positive, it becomes the
exponent of an exponential function. Thus muscle activation is modeled as
an exponential function of the difference between dynamic lambda and the
actual joint angle. This computation is carried out separately for the
flexor and extensor muscles. The "actual joint angle" as well as some of the
difference between the static and dynamic values of lambda presumable arises
from the output of the muscle spindles, sensors in the muscle sensitive to
muscle length and rate of change in length. The difference between the
dynamic activation threshold (represented as a joint angle) and the actual
joint angle is the error that is corrected in the model via negative
feedback - muscles are activated to produce torques in a direction that
reduces the difference.

BA: Muscle activation, E, is multiplied by a constant of proportionality to
produce the torque generated on the joint, which acts in opposing directions
for the flexor and extensor muscles. The difference in torques generates an
angular acceleration of the joint. This acceleration depends on the net
torque, divided by I (that's "i" not "one"), the moment of inertia of the
forearm. (This is the rotational analog to Newton's a = f/m (i.e.,
acceleration equals force divided by mass).

BA: Joint rotation accelerates under this net torque, this generates a
change in velocity at each time-step of the simulation dt, producing a new
velocity, and the change in velocity produces a new position of the joint.
The joint rotates until the torques balance. In the absence of an external
disturbance or encountering the limits of joint motion, the forearm moves to
the angular position specified by R and maintains a level of co-contraction
determined by C.

BA: I hope my explanation is clear enough without presenting the actual
formulae. These are available in the Lan and Zhu (2007) paper. The values
of all constants in the simulation were also taken from that paper, page
200.

Bruce

--
Richard S. Marken PhD
www.mindreadings.com

The only thing that will redeem mankind is cooperation.
                                                   -- Bertrand Russell