|Home | About | Journals | Submit | Contact Us | Français|
Properties of neural controllers for closed-loop sensorimotor behavior can be inferred with system identification. Under the standard paradigm, the closed-loop system is perturbed (input), measurements are taken (output), and the relationship between input and output reveals features of the system under study. Here we show that under common assumptions made about such systems (e.g. the system implements optimal control with a penalty on mechanical, but not sensory, states) important aspects of the neural controller (its zeros mask the modes of the sensors) remain hidden from standard system identification techniques. Only by perturbing or measuring the closed-loop system “between” the sensor and the control can these features be exposed with closed-loop system identification methods; while uncommon, there exist noninvasive techniques such as galvanic vestibular stimulation that perturb between sensor and controller in this way.
How do nervous systems transform sensory signals into motor commands for control? Engineering analyses, such as closed-loop system identification (system ID), promise to help answer this question. We focus in this paper on sensorimotor stabilization behaviors because they are particularly amenable to system ID. In these behaviors, animals robustly, automatically, and repeatably modulate muscle commands to drive sensory signals to desired equilibria (or limit cycles). A fitted model of an animal’s dynamical closed-loop response to a moving sensory scene can reveal features of the sensorimotor processing necessary to perform the stabilization behavior (Cowan and Fortune 2007).
Numerous systems in biology involve sensory guided stabilization of a mechanical system. For example, cockroaches follow along walls using their antenna by steering so as to stabilize the antennal tactile measurement (Camhi and Johnson 1999), which can be modeled as a sensorimotor stabilization task (Cowan et al. 2006; Lee et al. 2008). Weakly electric knifefish swim forward and backward to maintain their position relative to a longitudinally moving refuge (Rose and Canfield 1993) which can exploited to empirically determine the closed-loop sensorimotor transfer function (Cowan and Fortune 2007). Likewise, honeybees balance optic flow to remain in the center of a narrow passageway (Srinivasan et al. 1999, 1991; Si et al. 2003). Similar experimentally tractable sensorimotor stabilization behaviors occur in blowflies (Kalb et al. 2006; Balint and Dickinson 2004; Boeddeker and Egelhaaf 2005) and hawkmoths (Frye 2001; Sprayberry and Daniel 2007).
While our results are relevant to the aforementioned animal systems, the principle motivation and focus of our paper is bipedal standing balance in humans. Like the examples above, bipedal standing balance occurs in a closed sensorimotor loop (Horak and Macpherson 1996; Kuo 1995; van der Kooij et al. 2001; Kiemel et al. 2002; Carver et al. 2005) (see Fig. 1). For this system, we identify three components of this loop: the controller, N, the body, G, and the sensors, H. For human balance, N is the central nervous system, G is a multilink inverted pendulum (van der Kooij et al. 1999; Kuo 2005) actuated by muscles, and H is an array of vestibular, visual and somatosensory receptors (Horak and Macpherson 1996; Kuo 2005). The signals between pairs of “boxes” can be measured, such as with electromyogram recordings or force measurements, u0 (Fitzpatrick et al. 1996; van der Kooij et al. 2005), body sway, z0 (Oie et al. 2002; Kiemel et al. 2002), or neurophysiological recording, y0 (Ramcharitar et al. 2006). Of course, the available measurements will determine other modelling choices: for example if G includes muscle dynamics, then u1 cannot represent force. Likewise the signals can be perturbed, such as with muscle stimulation (Sponberg et al. 2008) or mechanical perturbation (Jindrich and Full 2002), wu, with visual scene motion wz (Oie et al. 2002; Ravaioli et al. 2005), or with a galvanic stimulus wy (Fitzpatrick et al. 1996). Of course, these perturbations must be delivered in closed loop, as the human or animal is performing the behavior under study.
By specifying known perturbations and using closed-loop system identification (Kiemel et al. 2002; van der Helm et al. 2002), researchers can create mathematical models for the individual components N, G, and H. Methods for identifying these components other than closed-loop system identification, for example “breaking the loop,” have been explored in other contexts such as the fly optomotor response (Heisenberg and Wolf 1988). These experiments have not been possible for human posture control because experiments with humans must obviously be noninvasive. Breaking the loop has several other disadvantages: in particular the subsystem under investigation may operate in a different (non-physiological) regime. Nevertheless, some related experimental manipulations with humans are used. For example, sway referencing uses feedback to attenuate ankle proprioception and/or vision (Nashner 1981; Nashner et al. 1982; Allum et al. 2002). But these manipulations fall short of opening the feedback loop because vestibular information remains present. Moreover, these experiments cannot be performed on patients with loss of vestibular function because the patients fall.
Posture control is often modeled as a linear system (Kuo 1995, 2005; Kiemel et al. 2002; van Soest and Rozendaal 2008). Define Wz, Wy, and Wu, to be the Laplace transform the inputs, respectively, wz, wy, and wu. Similarly, define Z0, Y0, and U0 to be the Laplace transform of the outputs, respectively, z0, y0, and u0. Then the transfer functions from input to output can be calculated and written as
where N = N(s), G = G(s), and H = H(s) are the transfer functions of the body, sensors, and controller, respectively.
We will show that under assumptions commonly made about posture control, zeros of N(s) cancel the poles of H(s). A common set of assumptions (stronger than needed) that imply these cancellations are (1) the system operates in the classical linear-quadratic-Gaussian context, (2) the nervous system implements optimal control based on a correct internal model of the plant, (3) the internal plant model includes a model of organism’s sensors, but (4) the optimality criterion does not penalize the sensor states.
In the classical linear-quadratic-Gaussian (LQG) context, a controller tries to regulate the setpoint of a linear system (described below using state matrices A, B, C and D) corrupted by Gaussian noise while minimizing a quadratic cost function (called J below). The optimal controller uses feedback from a optimal state estimator (a Kalman filter). The feedback gains (called K) depend on the optimality criterion but the optimal state estimator does not. The optimal state estimator depends on an internal model of the plant (with state matrices Â, , Ĉ, and ), as well as the Kalman gains L. The Kalman gains are determined by the internal model, and the statistics of the noise. A good reference with figures for optimal control in the context of human posture is Kuo (2005). A family of state estimators, called Luenberger observers, can be obtained by relaxing the requirement that the observer gains L, optimal‥ Our cancellation result holds for any choice of L. Moreover, our results do not require that the nervous system’s internal model be correct, but cancellations will occur only for the modes of the sensors that are correctly modelled in the nervous system. Our meaning of the phrase “correctly modelled” is subtle, especially in the case of multiple-input multiple-output systems, but it will be made precise below.
Under the scenario we have outlined above, in order to detect the poles of the sensor with closed-loop system identification, the system must be measured between H(s) and N(s) or perturbed between H(s) and N(s)—otherwise sensor dynamics are masked in the closed loop by the pole-zero cancellations resulting from the series connection N(s)H(s). More precisely, Proposition 2 shows that the pole-zero cancellations must occur in transfer functions for Z0 and U0 (i.e. in (1) and (3)). Thus, the only way to avoid the cancellation is to measure between H and N. Alternatively, however, perturbing between H and N may also prove useful: we show that in a simple example a double zero appears in the transfer functions from Wy to Z0 and from Wy to U0. Importantly, this double zero appears at the masked pole of H and one zero remains uncancelled in the closed loop. Thus, with a perturbation at wy sensory dynamics have an effect that is in principle measurable in the closed loop, even without measuring y0. These results suggest that to fully disentangle the respective contributions of sensory dynamics and downstream sensorimotor control, experiments must be designed to somehow access the signal between the two processes (Lee et al. 2008).
Suppose that body and sensor systems, G and H, are both linear, time-invariant, dynamical systems realized in state space as
Here, xb is the body state vector and xs is the sensor state vector; these vectors are of arbitrary (finite) dimension. We have dropped the subscripts on the signals u, z and y, because for the purposes of the present discussion (i.e. for deriving the controller with a zero that cancels a mode of H) we can assume that no perturbations are delivered: wu = wz = wy = 0, so that u0 = u1 = u, etc. The plant, the series connection between G and H, is given by the system
where is the concatenated vector of state variables. A simple calculation shows that
We adopt the increasingly popular view (Todorov and Jordan 2002; Bhushan and Shadmehr 1999; Kuo 2005, 1995; Kiemel et al. 2002; van der Kooij et al. 2001) that the animal’s neural controller N takes the form of (optimal) feedback from a state estimator. In particular, we presume a linear, time-invariant, observer-based neural control system:
where L is the observer error feedback gain, K is the state feedback gain, and (Â, , Ĉ, ) is nervous system’s model of the plant, H G in (6). Up to now, we have made no assumptions about the number of inputs and outputs of the various subsystems (except that they are internally compatible), the structure of the plant model, nor the optimality criteria used to select L and K.
The purpose of the remainder of this section is to illustrate, via a simple example, the more general result presented in Sect. 3. Toward that end, further suppose that the plant (6) is single-input–single-output (SISO) and that the actual sensor system H in (5) has a pole (perhaps one of many) at the value −α. The SISO assumption forces the signals u and y to have dimension one. We make no assumptions about the dimension of z, the signal between the body and the sensors.
For illustrative purposes, assume the observer is based on the following model of the plant:
where , , are all scalars, i.e. the observer makes the possibly erroneous assumption that the plant is second order. We notate the state space matrices for the observer’s model of the plant with hats to distinguish them from the actual state space matrices of the plant—we allow the two sets of matrices to be different. The components of the state of this internal model are the estimated body state b and estimated sensor state s. We allow the estimator’s model to be incorrect both in parameters and in number of states. In this SISO example, the internal models of both G and H are first order, but we emphasize that we have made no assumptions about the orders of the true G and H—they could potentially be much larger than one. The zero that appears in the upper right corner of Â follows from the series connection between G and H. Our choice of internal model implies that the controller knows the plant has this structure. Finally our internal model of the plant assumes that the sensor system has a pole at the value −.
In this illustrative example, the observer-based control (9) reduces to
We do not yet specify the elements of the matrices K and L, however compatibility with our previous assumptions determines their dimensions:
We make two additional assumptions, both also having analogues in the general result. First, we assume ks = 0, or in general, that the feedback does not depend upon the estimated sensor state variable. This assumption, commonly made by postural modellers (Kuo 2005; van der Kooij et al. 2001) will be justified in the next section (Proposition 1) as a consequence of optimal control with a cost function that involves only mechanical states and not sensor states. Second, we assume = α, namely that the sensor system H actually has a pole where the internal model assumes.
We now verify that the open-loop transfer function of this controller is given by
for some constants η1 and η2 that depend upon all parameters of the controller (including ).
To prove this statement, we write the matrix F in terms of its components:
Now we use this representation of F, and the formula for the inverse of a 2 by 2 matrix, to expand (14):
Using that ks = 0, we expand the second column of F:
It can now be verified that that the numerator of N(s) is kblb(s + ) and that the denominator of N(s) is s2 − s(F11 + F22) + (F11 F22 − F12 F21), as claimed.
This transfer function has a single zero at −. But we have assumed that H has a pole at −α and that the internal model knows this, i.e. = α. A zero in N coinciding with a pole in H implies a pole-zero cancellation in the closed loop when the measurement is taken at u0 or z0. Thus, the controller masks the sensor pole at −α. Our masking result follows from a few simple assumptions that allow the internal model to be wrong in many respects.
Even if the signal cannot be measured between H and N, there is some hope of identifying the sensor dynamics from perturbations between these two blocks. Consider the transfer functions from Wy to U0 and Wy to Z0, (embedded in (3) and (1), respectively). Write
Because H(s) and N(s) are rational transfer functions and each has a factor of (s + α) (in, respectively, the denominator and numerator) (s) and Ñ(s) represent, respectively, H(s) and N(s) separated from the factor (s + α). The purpose of this manipulation is to demonstrate the cancellation between the factors of (s + α) and, in particular, where those factors remain uncancelled. With these substitutions and setting wz and wu to zero yields the following transfer functions
Thus, an uncancelled zero appears at the (open-loop) sensor pole in the transfer functions from Wy to Z0 and from Wy to U0.
We assume G(s), H(s) and N(s) are given by (4) through (9), except that now we make no assumptions about the dimensions of u, y, xb, xs or . In particular, we allow the plant and internal plant model to both be multiple input and multiple output (MIMO). The internal plant model can still be wrong both in parametrization and in numbers of states, but we assume concatenates estimated states for the body, b, with estimated states for the sensors, s.
Recall that the plant refers to the series connection of the body, G, and sensors, H, which has the following block structure:
The nervous system’s internal model of the plant has an analogous block structure implied by the series connection between G and H, known to the nervous system. Specifically,
We reiterate that neither the elements, nor the dimensions, of the internal plant model’s matrices (24) and (25) need be the same as the corresponding quantities for the actual plant model, (22) and (23), except that the upper right zero block of the plant A matrix is known by the nervous system.
As before, we assume N(s) implements state feedback based on a Luenberger observer (9) with the form
As before, F = Â − ( − L ) K − L Ĉ. Moreover we assume that the dimensions of L and K are compatible in (26). These matrices have block structure:
The matrix L can chosen so that the observer is a Kalman filter, making estimate the state x optimally with respect to the internal model and the assumed statistics of Gaussian noise. Moreover, the matrix K can be chosen to provide optimal control with respect a quadratic criterion, making the observer an optimal linear-quadratic-Gaussian controller. We do not explicitly represent noise in our equations, but we do allow for noise. Specifically, noise could corrupt the signals wy, wu, and wz, in Fig. 1. If the internal model is correct and K and L are chosen appropriately (with respect to a quadratic criterion), then in the presense of the noise, our controller performs better than any other with respect to the chosen criterion.
We emphasize that our results do not depend on the assumption that L and K are chosen optimally. However, as in the SISO case of Sect. 2, we do assume that there is no feedback directly from the sensory states, namely Ks = 0. This would be the case assuming a quadratic optimal controller with no penalty on the sensory states, as established by the following result.
Suppose K is the unique feedback gain that optimizes the quadratic cost function
where Q and R are symmetric matricies, Q is positive-semidefinite, and R is positive-definite. Suppose further that J does not depend upon sensor variables, (A, B) is stabilizable and (A, Q) is detectable. Then Ks = 0.
The feedback gain K satisfies
where S is the unique symmetric positive-semidefinite solution to the algebraic Riccati equation (Bryson and Ho 1975):
The block structure of x (body, sensor) implies a block structure for Q. Because the objective function does not depend upon the sensor variables, Q has the following block structure:
Note that this assumption would follow from a hypothesis that postural control system has evolved to stabilize the body, not the sensors. Given this structure for Q, a simple calculation verifies that S is given in block structure by
where Sb is the unique positive semidefinite that satisfies
and where (AG, BG, CG, DG) are the system matrices for the body G(s). The optimal feedback gain is then given by [Kb 0] where
Finally, note that the result Ks = 0 still follows if the system matrices used in this derivation from the internal model of the plant, rather than the actual plant. This generality holds because we have assumed that Â has the block structure determined by the series connection of G and H.
In the SISO case, if − is a pole of the observer’s model of the sensor (i.e. an eigenvalue of Âs) then the controller N(s) has a zero at −. Now consider the MIMO system and suppose − is an eigenvalue (real or complex) of Âs. Under what conditions does N(s) have a transmission zero at −?
We refer the reader to Levine (1996) for the following definition: we say that s is a right zero of the system (with n states, m inputs, and p outputs)
if vectors ξk and uk exist, both not zero, so that
The matrix on the left hand side of (37) is called the Rosenbrock system matrix (RSM).
We have assumed that there exists a vector such that Âs = −. Let q = Ĉs. Direct calculation shows that the vector
lies in the null space of the RSM
if and only if s = −. This result shows that the controller N(s) has a transmission zero at .
Now we show that a pole-zero cancellation occurs in the closed loop:
Let be the set of available perturbation inputs to the plant and let yf = [zT, uT]T (with no direct measurement of y) be the set of available measurement outputs. If
then there is a pole-zero cancellation, i.e. α is a pole of the closed-loop system but the mode corresponding to α is unobservable.
In less precise language, we can rephrase Item 4 as the condition that the internal model’s observation of the sensor mode is correct. Item 4 will always be satisfied if the system has one observation (i.e. if C and Ĉ have one row), and also if the quantities Csυ and Ĉs are nonzero. Alternatively, Item 4 will always be satisfied if the internal model of the sensors is correct. Finally for the closed-loop system to remain stable we must also assume that the canceled modes of the sensors have negative real parts.
Let xf = [xT, T]T be the full system state. Then the full system can be written as
the block structure of Bf is unneeded for the proof, and
Direct calculation verifies that the vector
is an eigenvector of Af with eigenvalue −α. (Here 0 indicates a zero vector of appropriate dimension.) Thus, the sensor pole is a pole of the closed-loop system. Nevertheless Cfυf = 0 indicating that the mode associated with α is unobservable due to a pole-zero cancellation.
We have shown that pole-zero cancellations between the sensors and the controller occur under modeling hypotheses commonly made about the human postural control system. For examples, our results hold under the classical optimal linear-quadratic-Gaussian paradigm where a (correct) internal model of the plant includes the sensors, but the cost function for optimality does not penalize sensor states. For clarity, we concede that our results do not prove rigorously that these pole-zero cancellations also appear in biological systems. Indeed, we cannot establish that all the assumptions needed for the proof hold in nature. For example, the subsystems, H, N and G, are unlikely to respond linearly to all possible signals. Instead, our results establish that mode-masking, while it may naively seem pathologically unlikely, actually appears generically in sensorimotor models under hypotheses considered reasonable. Our results suggest that the possibility of such masking in biological systems should not be assumed to be negligible. Thus, to reveal the dynamics of the sensorimotor chain, methods of identifying the subsystems should be designed accordingly. To unmask hidden dynamics these methods must include measuring, perturbing or breaking the loop between the subsystems.
A special case of sensory dynamics is a time delay. Time delays exist in all sensor systems, biological or otherwise. A linear sensor in series with a time delay is a linear dynamical system, however it is not one that can be placed into the finite dimensional state space format, (22) and (23), that we have assumed to prove our results. Specifically, a time delay is an infinite dimensional dynamical system.
To make this discussion concrete, suppose that the sensor system H behaves as a first order linear system with a pole at −α followed in series by a time delay. Now, consider carefully the proofs about the appearance of the zero at in the SISO controller or transmission zero at in MIMO controller. Remember that is the value of α assumed by the nervous system. Notice that the only assumptions used in these proofs concerned the internal model of the plant and not the actual plant. If the internal model of the plant ignored the time delay, but satisfied our other hypotheses, it would still have a (transmission) zero at −, which could cancel the mode at −α. For this cancellation to occur we must of course have match between sensor pole and its internal model: −α = −. Also, for MIMO systems, we must have a match in the observation of this mode (Item 4 in Proposition 2). But the presence of infinite dimensional dynamics in the plant does not prevent it. Thus the mode masking phenomena we study here can occur even in the presence of time delays. But only the finite number of modes that are assumed and modelled correctly by the nervous system are cancelled by the controller—the rest remain.
On the other hand, the internal model of the plant might not ignore the time delay. One way the internal model might account for the delay is with a finite dimensional approximation. Such an approximation can account for delay arbitrarily well if the order is sufficiently high. In this situation, assuming the other hypotheses are satisfied, the controller has (transmission) zeros at the poles of its internal delay approximant. Thus the controller will mask a finite number of the modes of the delay. However this does not mean that the controller is masking the delay, even approximately, because the controller does not invert the transfer function of the delayed sensor. In particular all of the zeros of the delayed sensor remain in the closed loop.
We raise the possibility that the mode-masking phenomenon that we study here, could, in the future, be supported experimentally. We discuss this possibility in the context of posture control and the vestibular sense. The primary afferents of the semicircular canals have been modeled as high-pass filters of velocity with a time constant of about 5 seconds (Mergner 2002; Wilson and Melvill Jones 1979). This time constant arises from the biophysics of the end organ. More detailed models place an additional time constant significantly faster (0.003 s, for both the regular and irregular units) as well as one much slower (32.4 s, for the irregular units) (Wilson and Melvill Jones 1979).
In engineering practice the dynamics of sensors often occur in a higher frequency regime than the frequency regime of the dynamics of interest. In this scenario, an engineer might not penalize the sensor states in the cost function for optimal control, just as we have assumed. However contrary to our assumptions, an engineer might not even include the sensors in the forward model (called the “internal model,” above), because they are not needed for control, and contrary to our results, the sensor’s high frequency modes would remain uncancelled in the closed loop. But the time constants of the vestibular sensor are slow compared to engineered sensors and human postural sway has a prominent time constant of roughly the same magnitude (5–10 s) as the 5 s time constant of the vestibular afferents. This time constant persists in the postural sway of patients with profound bilateral vestibular loss (unpublished observation), suggesting that its origin is not vestibular in nature. The omission of the vestibular state in the forward model would severely impair the hypothetical engineer’s controller.
Because of the frequency overlap, and because the vestibular time constants are known fairly precisely, we propose the 5 s vestibular time constant as a good candidate for experimentally verifying our mode-masking predictions. Specifically, we propose searching for a coinciding zero in the neural controller. One way to detect such a zero, without breaking the loop, is to use closed-loop system identification with a perturbation that affects the vestibular signal after the dynamics of the sensor. There is an experimental technique which accomplishes such a perturbation: galvanic vestibular stimulation (Fitzpatrick et al. 1996). The challenge of using such a stimulus with system identification is controlling the stimulus precisely. A wider range of tests may be possible with animal models, such as postural balancing in cats (Ting and Macpherson 2005; Lockhart and Ting 2007). As sensorimotor system ID techniques mature, we expect that the questions raised in this paper will be addressed in a variety of animal and human experiments.
We wish to thank Jim Freudenberg for carefully reading the manuscript and providing helpful comments.
This work was supported by the National Science Foundation under Grant No. 0543985, and by the National Institutes of Health under Grant No. 2RO1 NS35070.
Sean G. Carver, Department of Psychological and Brain Sciences, The Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218, USA.
Tim Kiemel, Department of Kinesiology, School of Public Health Building, The University of Maryland, College Park, MD 20742, USA.
Noah J. Cowan, Department of Mechanical Engineering, The Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218, USA.
John J. Jeka, Department of Kinesiology, School of Public Health Building, The University of Maryland, College Park, MD 20742, USA.