The maximum isometric force *Fmax* for elbow muscles was constrained to be smaller than that of the shoulder muscles, since the latter need to support and transport a larger mass than the former. Also, in the literature the strongest shoulder muscles (such as the deltoid) are consistently reported to be stronger than elbow muscles (Nijhof and Kouwenhoven, 2000; Garner and Pandy, 2003; Holzbaur et al., 2005). The strength difference in our optimized model is rather big. For the purpose of our investigation, however, this would pose a problem only if optimized muscles were unrealistically strong, and if neural control would exploit this strength to overpower the interaction torques at the elbow joint. But this is not the case. As we demonstrate below, the combination of muscle strength and activation levels leads to a realistic range of dynamic torques throughout the movement. For example, interaction torques can reach the same level as muscle torques, and maximum shoulder torques (on the order of 10 Nm), are a multiple of maximum elbow torques (about 2 Nm), similar to measures from human subjects (e.g., Galloway and Koshland, 2002).

Finally, maximum contraction velocities (*vmax* approx. 11.4 *l*0/*s*) fall within a physiologically plausible range. Zajac (1989), for example, assumes an average of about 10 *l*0/*s*; Ranatunga (1984) measured values between 7 and 13 *l*0/*s*; for the empirically based model used in Kistemaker et al. (2007) *vmax* is not specified numerically, but visual inspection indicates values about or greater than 10.

The described morphology is only one of a wide range discovered by the search procedure. Others were found with similar performance but great variation in selected parameter values, indicating that the task underspecifies the required musculoskeletal properties.

#### **3.2. MOVEMENT KINEMATICS**

To assess whether the optimized model reproduces empirically observed kinematic invariants (Morasso, 1981; Soechting and Lacquaniti, 1981; Atkeson and Hollerbach, 1985; Flash, 1987; Flanagan et al., 1993), in **Figure 4** we show movements performed by the best optimized spinal circuit. Trajectories are approximately straight with slight curvature, feature small hooks at the target positions, and exhibit the characteristic bell-shaped velocity profile. This is true for all four movements, i.e., independent of the direction of interaction torques or starting posture.

In panel B, we plot example joint trajectories for movements *WA* and *RA*. Consistent with earlier findings (Ghafouri and Feldman, 2001), we observe that the commanded trajectory (dashed) is significantly shorter than the movement period. The optimized duration covers 44.6% of the actual movement. Also note that the commanded (and desired) trajectories of individual joints are offset slightly in time, a result of their derivation through inverse kinematics from a hand path planned in extrinsic space.

#### **3.3. MOVEMENT DYNAMICS**

The demonstrated kinematic invariants alone are not sufficient to imply active compensation, or accommodation, of interaction torques by the spinal cord. Although we know that central commands in our model carry no anticipative corrective components—the threshold shift is always of the same monotonic form—we need to show that the spinal cord can transform these identical control signals into muscle activation patterns that differ qualitatively with the direction of interaction torques (Cooke and Virji-Babul, 1995; Latash et al., 1995; Gribble and Ostry, 1999; Galloway and Koshland, 2002; Debicki and Gribble, 2005). **Figure 5** shows torque patterns produced by our model for movements starting from initial posture *SA* (corresponding to kinematics shown in **Figure 4B**).

Firstly, we observe that the interaction torque experienced in one joint (dashed) strongly correlates with the total torque (filled) in the other. Secondly, comparing whipping and reaching movements, we find that elbow torques vary with the direction of shoulder motion even though elbow kinematics are held constant. In whipping movements interaction torques due to shoulder movement initially oppose the movement of the elbow. This effect is significant, as peak interaction torque in the elbow is slightly greater than the muscle torque. Interestingly, in this case the two torque components almost cancel out, which leads to a delay in the onset of elbow motion of exactly the duration necessary to follow the desired hand trajectory. Whipping movement *WB* (not shown) differs in this respect, in that the interaction torque here is slightly smaller than muscle torque, resulting in no such onset delay (again, as the desired hand movement requires). The shoulder, in contrast, is subject to only minimal interaction torques, and its movement is consequently dominated by muscle torques. This was also found to be the case for the second starting posture.

During reaching movements, in comparison, our model shows that interaction torques at the elbow initially assist the motion. They are equal in sign to the muscle torques and on the same order of magnitude. Though this effect is stronger in the elbow, interaction torques also assist shoulder motion, but in this case are significantly smaller than muscle torques. We thus find that both types of movement are characterized by a shoulder-centered pattern, in which initial shoulder motion is generally dominated by muscle torques, while elbow motion is produced with significant contribution from interaction torques.

The pattern of opposing and assisting effects of interaction torques also determines the overall muscle effort required. In the opposing case, elbow muscle torques need to be larger than in the assisting case to compensate for interaction torques, even though the kinematics of the motion is essentially the same (a 30◦ flexion) in both cases. This is particularly salient in elbow torques during reaching. Here the muscles do not contribute at all to the braking forces that terminate the movement, i.e., we observe no forces resulting from antagonist activity. The braking pulse in this case is exclusively due to interaction torques created by shoulder motion.

In short, the movement kinetics in our model confirm that the spinal circuitry successfully transforms simple descending control signals into muscle force patterns that qualitatively differ with the direction of interaction torques in such a manner as to accommodate them.

#### **3.4. NEURAL DYNAMICS**

In **Figure 6** we plot the activity of α-MNs and the individual muscle torques they provoke (taking into account muscle activation dynamics and changing moment arms) for the same movements as shown in **Figures 4**, **5** (*WA* and *RA*). Neural activity exhibits the characteristic bi- or tri-phasic burst patterns observed empirically (Ghez and Gordon, 1987), i.e., we can generally identify an accelerating agonist burst followed by a decelerating antagonist burst, and sometimes a third agonist burst arresting the motion.

Not surprisingly, given the torque patterns described above, elbow muscle activity varies with the direction of motion in the shoulder, despite virtually identical elbow kinematics. Muscle activity is generally greater when both joints move in the same direction, i.e., when interaction torques oppose the movement (whipping). For example, integrating over the corresponding area under the curve, we find that motor neuron activity associated with the first elbow extensor burst (which is the agonist in these movements) is approximately 0.033 for the whipping movement, but three times smaller (0.011) for the reaching movement. The same is true for the antagonist (here the elbow flexor). In the former case the area of its burst is 0.024, and when interaction torques assist the motion it vanishes completely.

In contrast with some empirical data, we observe in our model no systematic time lag between onsets of activity in agonists acting on different joints; in particular, we find no temporal organization from proximal to distal joints (Karst and Hasan, 1991; Gribble and Ostry, 1999). Even though such a lag seems to be present in movement *RA* (in **Figure 6** compare first α-MN bursts of elbow and shoulder in the reaching condition), this was not the case for all reaching movements.

For reasons of space we do not present a full analysis of the neural dynamics exhibited by the optimized spinal circuitry. We suggest, however, an explanation for the suppression of the elbow antagonist burst during reaching movements. In the optimized spinal circuit, the antagonist burst (and its suppression) can be traced back to two opposing influences on its α-MN pool. First, spindle feedback excites the α-MN in proportion to deviations from desired position and velocity. But secondly, spindle activity also drives IaIn activity, and via connections from IaIn to homonymous IbIn (Jankowska, 1992) also the latter interneurons, which inhibit α-MN activity. A second inhibiting influence in the optimized network originates in IbIn intersegmental connections from muscles acting on the other joint. The size and shape of the antagonist burst therefore is the result of a balance between spindle feedback and IbIn activity, which in turn is modulated by Ia interneurons. In whipping movements, presumably because interaction torques initially oppose the movement, position and velocity errors initially grow relatively large, resulting

in spindle feedback sufficient to overcome the aforementioned inhibiting factors. In reaching movements, on the other hand, interaction torques partially "do the work," which leads to smaller deviations from the desired state, and hence spindle feedback that is more easily suppressed by the same inhibiting factors.

almost exclusively due to interactions torques produced by the shoulder.

#### **3.5. GENERALIZATION**

The model presented above has been optimized for movements of a certain amplitude and speed (in joint space) and in two separate regions of the arm's workspace. Here we briefly address whether it generalizes to other types of movement that it has not been optimized for. We start by testing the model's performance as we shift the two starting postures in joint space by values ranging from −25◦ to +10◦ (elbow and shoulder angles are shifted by the same amount), while keeping duration and amplitude fixed at the original values. As **Figure 7A**

**FIGURE 7 | Model generalization ability.** Performance of the optimized controller as the two starting postures are offset by a given number of degrees in joint space **(A)**; as movements vary in duration **(B)** or amplitude **(C)**; and as duration and amplitude are changed in proportion such as to maintain average velocity **(D)**. Measures on horizontal axes are relative to values used for optimization, which are indicated by vertical lines. Performance drops the more movements differ from the optimized kinematics. In red we plot performance when cocontraction signals are adapted separately for each desired movement.

demonstrates, the optimized controller shows specificity for the area of the workspace encountered during optimization. Performance decreases in both directions from the original postures.

Similarly, performance drops quickly as we change the duration of the movements to be shorter or longer than the original (panel B), or as the amplitudes are increased or decreased relative to those used during optimization (panel C). If we scale amplitude and duration equally, such that the average velocity remains constant, the drop in performance is less dramatic, but the overall picture is the same (panel D).

In the above tests, none of the system's control signals were re-adjusted for the varying movement kinematics. It cannot reasonably be expected, however, that the resulting movements should be well executed if, for example, the amount of cocontraction (the component λ*co*) is not tuned to the speed demands of the desired movement. It has been shown in human subjects, for example, that movement velocity correlates with muscle cocontraction (Gribble et al., 1998). Also, the balance of open-loop antagonist muscle activity specifies the equilibrium of the system in statics, implying that a wrong selection of this balance (given a specific target) could lead to spinal circuits and musculotendon system driving toward incompatible equilibria. At least these components of the feed-forward command therefore have to be chosen selectively for each particular movement. To test whether this is indeed sufficient to achieve reasonable performance, we first re-optimized all parameters of the original model after adding two new movements to the performance evaluation: the first is identical to *WA*, except that its duration is 0.4 s instead of 0.3 s; the second one is also similar to *WA*, but here both the duration and amplitude are 20% larger, such that the average velocity remains the same. The purpose of these additional evaluations is to avoid optimization of controllers that are overly specialized on the speed and amplitude demands of the four original movements. In a next step we then choose a few movements that the original system performs badly on and re-optimize only the cocontraction signals.

In **Figure 7D** we compare performance over a range of movement amplitudes (but constant average velocity) when the cocontraction signals are re-optimized (red) with those not adapted for each specific movement (black). In the former case performance remains almost constant, dropping no more than approximately 2% (from 98.9 to 97%). The performance of the non-adapted system, in contrast, drops to only 54%. We also tested a few of the other conditions presented in **Figure 7**. For example, the performance when desired joint rotations are reduced by 10% (while keeping the original duration fixed, thus leading to reduced average velocity), is 98.3%, compared to only 8% in the non-adapted system (see panel C). Equally, if the desired duration of the movement is increased by 20% while keeping the amplitude the same, performance of the adapted system is 98.7%, compared to only 67% for the non-adapted system. Performance is improved further if in addition to cocontraction signals we also tune spindle sensitivities to the desired movement, i.e., when we adapt the properties of the gamma pathways such that the strength of position and velocity feedback depends on the desired amplitude or velocity of the movement (data not shown).

To test the model's capacity for producing movements not only of different amplitudes and durations, but also in different directions, in **Figure 8** we show the performance of the spinal circuit when optimized for an increasing number of movement directions. All movements here have a duration of 0.3 s and follow desired center-out hand paths 10 cm in length. Generally, hand paths are essentially straight (panel A, MED averaged over four movements: 1.8 mm). When the number of directions increases (panels B and C), some paths remain almost perfectly so, while others begin to show slight curvature (average MED for six movements: 2.3 mm; for eight movements: 4.2 mm). Almost all paths

movements have an amplitude of 10 cm, a duration of 0.3 s and proceed from the center (black dot) to targets spaced equally along a circle. Hand paths for most movements are approximately straight, except for those in the directions of 90◦ and 270◦ (dashed), which reflect a limitation of the biomechanical model.

resemble the variation of curvatures observed in human reaching movements, except for two movements in panel C (toward targets located in the directions of 90◦ and 270◦). Their large curvature and inaccurate termination reflect biomechanical constraints due to our muscle model, as these movements could not be optimized successfully even in isolation (average MED without these movements: 1.8 mm). This is indicative of the limitations of our simplified musculoskeletal system.

We also note that all of the optimizations described in this section were performed without the velocity error term in the threshold model, and without the power transformation of the viscosity term. Further tests, which we do no present here in detail, show that a control model without spinal interneurons depends on the presence of at least one of these terms to approach the performance of the spinal circuit.

In summary, when control parameters are tuned to the kinematics of the desired movement (e.g., cocontraction is matched to desired velocity), smooth movements can be performed accommodating interaction torques. This is the case for movements of different amplitude and velocity, in different areas of the workspace, and in different directions.

#### **3.6. THE ROLE OF IBIN ACTIVITY**

As mentioned in the introduction, there is reason to speculate that intersegmental force feedback may provide a mechanism by which the nervous system compensates for interaction torques during multijoint movements. One can in fact conceive of three different roles: force-related feedback could modulate descending commands by acting in long loops through the CNS; coordinate muscles acting on different joints through intersegmental neural connections; or act indirectly through the mechanical coupling of joints and the sensed effect of internal loads in a manner akin to non-neural coordination mechanisms in the stick insect (Schmitz and Stein, 2000; Cruse et al., 2007).

We investigate here the latter two options by performing simulated ablation experiments in which we re-optimize the system's parameters after removal of different parts of the spinal circuit. **Figure 9** shows hand trajectories for the worst of the four movements in each ablation condition, which happens to be *WA* in all cases. Optimizing the system without intrasegmental Ib connections reduces the system's performance to a relatively small degree (compare panel A and B). The curvature of the trajectory becomes more pronounced and the movement is arrested less effectively. Movement error, measured by MED, increases from approximately 3.5 to 6.6 mm (and from 2.37 to 3.14 mm when averaged over all four movements). When Ib interneurons are completely removed from the network (panel C), curvature and endpoint oscillations become even more salient (MED = 12.9 mm; average MED = 6.24 mm). In addition, the time course of the trajectory deviates significantly from the reference. It initially leads the desired movement, but then fails to reach the required velocity and falls behind toward the end. For comparison we show in panel D the result when all interneurons are removed from the network, leaving only the threshold model, i.e., α-MNs and direct proprioceptive feedback. Although the trajectory in space shows no more curvature than seen in the other conditions, its time profile deviates much more

strongly from the reference, which is most salient in the velocity profile.

#### **3.7. MODEL SENSITIVITY**

To assess the sensitivity of the model we evaluate its performance over a range of deviations from the solution optimized for the four original movements (*WA*,*B* and *RA*,*B*). The optimization procedure encodes all parameters in a vector with component values in the range [0, 1]. For a given level of deviation μ we add random perturbations chosen uniformly from the range [−μ, μ] to all components, i.e., we select a random vector from within a hypercube of size μ centered on the optimized parameter vector. For every level of deviation we sample 100 such vectors and measure their mean performance. The result is shown in **Figure 10**.

It can be seen that performance drops gradually as the amount of deviation increases, indicating that the procedure has not found a "needle in a haystack." The relative smoothness of the error surface suggests that other optimization methods, for example those based on trial-and-error or gradient descent, should also be able to find good solutions. Gradient-based learning using trajectory errors has in fact been demonstrated in a similar model of movement control using spinal-like neural networks (Raphael et al., 2010; Tsianos et al., 2011). On the other hand, performance starts decreasing immediately. The absence of a plateau around the optimal set of parameters suggests that there is not a great variety of different models leading to the same performance. However, since we are only probing (an increasing) neighborhood of one optimal solution, we cannot rule out that other solutions with similar performance exist in other regions of parameter space (in which case local optima might in fact hinder the performance of gradient-based algorithms).

# **4. DISCUSSION**

The dominating view in motor control today suggests that the CNS controls the body using intricate internal models of its kinematics and dynamics in order to predict and directly control the muscle forces required to perform a desired movement. An alternative view, expressed in the equilibrium point hypothesis, proposes that the combined dynamics of spinal circuitry and musculoskeletal system provide a level of abstraction in the control hierarchy that allows the CNS to plan and control movements without requiring a representation of complex bodily dynamics. In this paper we provide evidence that this is plausible. Instead of anticipating upcoming interaction torques and adjusting central control signals accordingly, our model suggests that the CNS may in some cases—such as the type of reaching movements considered here—be ignorant of musculoskeletal dynamics and offload the coordination of muscle forces required for a particular, kinematically defined, movement goal to circuitry at the spinal level. We do not rule out that prediction or implicit anticipatory mechanisms might be involved in other cases, such as faster or more complex movements.

Several studies (e.g., Almeida et al., 1995; Gottlieb et al., 1996; Gribble and Ostry, 1999) have found that during human limb movements EMG activity in muscles acting on one joint correlates with interaction torques arising from motion in another, and is often timed such that it precedes the onset of movement. These findings have been taken to imply that central motor commands are adjusted predictively to compensate for interaction torques. But it need not be true that any such adjustment takes place on the level of central control signals, nor that any form of prediction is involved. Firstly, it is clear that EMG activity has to vary systematically with upcoming interaction torques. If it were not to we would not observe hand paths that are approximately straight. Also, some muscle activity has to precede the movement, as it is necessary to initiate it. Furthermore, since we do not yet have sufficient knowledge about the precise nature of descending control signals, which only after integration with afferent and interneuronal signals results in observed EMG activity, it is impossible to conclusively deduce from current empirical data whether central commands are already adjusted for interaction torques, or whether they are transformed at a lower level to this effect. The model presented here demonstrates that accommodation of intersegmental loads on the spinal level is possible. Also, the fact that the temporal order of muscle activity—including prior to movement onset—seems to be relatively fixed and organized such that agonists in proximal joints precede distal ones, may reflect a more general organization of the control hierarchy, rather than specific and detailed predictions of upcoming dynamics.

The muscle and reflex model developed here involves significant simplifications when compared with the problem of controlling multijoint arm movements in the human body. For example, it does not take into account the effect of biarticular muscles, tendons or gravity. Spinal interneurons are modeled as simple leaky integrators, and fusimotor drives are represented only implicitly through the λ-model. Also, the Hill-type muscle model is employed here in its simplest form, ignoring, for example, the effect of calcium sensitivity on the force-length characteristic (Kistemaker et al., 2007), or the dependence of maximum shortening velocity on activation level (Chow and Darling, 1999). While the inclusion of tendons may be important, for example, when studying the physiological mechanism underlying the estimation of joint position (Kistemaker et al., 2012), we have argued in Section 2.2 that its effect on upper arm dynamics is limited. And if biarticular muscles are particularly well suited for internal load compensation (Gritsenko et al., 2011), then their inclusion in a model such as the one presented here should be expected to make the problem of feedback compensation easier. We also note that the problem of interaction torques in multi-segment limbs is in principal independent of how the joints are actuated, as long as the mechanism of actuation is not infinitely stiff (the problem hence also arises, for example, in any type of compliantly actuated robots). We therefore believe that none of the simplifications introduced in the model interfere with the basic goal of this study, which is to demonstrate that a single spinal-like neural network can transform simple descending control signals into muscle activation patterns that differ qualitatively with the direction and magnitude of interaction torques in a manner that is appropriate for the generation of smooth and straight hand trajectories. Moreover, the model achieves this with muscle and interaction torque patterns comparable to those observed empirically, and without the need for inverse dynamics calculations, prediction of upcoming loads, or having to learn adjustments to control signals for each individual movement. Nevertheless, we believe that further work aimed at testing the proposed mechanism in more realistic models is required to show how and whether the control scheme illustrated here is in fact employed in the case of human reaching movements and to generate detailed predictions that can be tested empirically and that go beyond the demonstration of a functional role for proprioception in intersegmental coordination (Ghez and Sainburg, 1995; Sainburg et al., 1995).

Despite the simplifications just mentioned, the model presents a significant complexity. This is mostly due to the need to represent spinal circuits explicitly in order to do justice to the hypothesis being investigated. Though the modeled circuits are still far simpler than real spinal circuits, we do not claim to have found the simplest or most efficient model able of interaction torque compensation, which was not our objective. In any case, the work presented here also demonstrates that the chosen methodology is indeed practical for asking questions about the potential functionality of spinal circuits, despite their complexity.

In addition to demonstrating the feasibility of feedback compensation for interaction torques in general, our model reproduces several features of reaching movements performed by humans. As demonstrated in Section 3.2, the hand trajectories generated are approximately straight and exhibit bell-shaped velocity profiles (Morasso, 1981; Flash, 1987; Flanagan et al., 1993). This is not surprising, since the movements were optimized to match minimium jerk profiles (Flash and Hogan, 1985). Exceptions from this mathematical ideal, such as uni- or bimodal curvature and hooks at the endpoints, are also in correspondence with empirical data (Flash and Hogan, 1985).

Regarding muscle and interaction torque patterns, Gribble and Ostry (1999)report that muscle torque applied at the elbow varies with the direction of shoulder motion across movements in which elbow kinematics are held constant (and, conversely, shoulder muscle activity varies with the direction of elbow motion when shoulder kinematics are constant). These dependencies are found even when active movement in only one joint is required (but the other is free to move) or when one of the joints is fixed (Debicki and Gribble, 2005). The direction of this covariation depends on whether the two joints move in the same or opposite direction. When shoulder and elbow move in the same direction, interaction torques arise in each joint that initially oppose that joint's intended motion and muscle torques in each joint increase with the emerging interaction torque. In contrast, when the joints move in opposite directions, interaction torques initially assist the desired movement, and muscle torques decrease proportionally. As we have shown in Section 3.3, our model behaves in the same way. Simulated elbow torques vary with shoulder kinematics, and overall effort depends on whether interaction torques are assistive or opposing. As reported in Section 3.4, this effect can be strong enough to completely suppress the normal antagonist burst in the elbow. In this case it is the interaction torque at the elbow alone that arrests the motion. Similar effects can be observed in human subjects, where in some cases assisting interaction torques reach levels such that an initial counteracting antagonist burst is required, instead of the more usual movement-initiating agonist activity (Cooke and Virji-Babul, 1995). Our model also exhibits a shoulder-centered pattern, where elbow dynamics are the result of approximately equal muscle and interaction torque contributions, while in the shoulder muscle torques are generally greater than the passively arising loads (at least initially). This is consistent with some empirical data (e.g., Galloway and Koshland 2002), although it has not been reported by Gribble and Ostry (1999).

With regard to the timing of muscle activity, we find no systematic temporal organization of agonist onsets from proximal to distal joints (Karst and Hasan, 1991; Gribble and Ostry, 1999). Investigation of a larger range of movements would be necessary to identify whether the model does or could exhibit such a strategy. We do, however, observe instances of onset delays in joint rotations (Karst and Hasan, 1991; Virji-Babul and Cooke, 1995), which result from the interaction of muscle forces and internal loads, and from the planning of movements in extrinsic space.

The optimized spinal circuits produce desired hand trajectories over a range of different amplitudes, speeds and directions (Section 3.5), with the exception of movements along a single direction, which was shown to be a limitation of the biomechanical model rather than its control. Further work would be required to determine whether the addition of biarticular muscles, for example, would broaden the operating range of the system, or whether the optimization procedure simply failed to identify a more capable biomechanical setup given the model's constraints. For the model to achieve good performance in all movement conditions, we have shown that some control parameters (such as coactivation level, or muscle spindle sensitivity) need to be selected on the basis of the desired movement speed and amplitude. But crucially, no details about the dynamics of the movement need to be known. The control signals are always simple—a monotonic shift in muscle reflex threshold and a constant level of coactivation throughout the movement—and do not depend on anticipated interaction torques or other aspects of the system's dynamics. The results suggest some further questions, however. For example, what are the *minimal* changes in control signals that allow for the control of movement speed or amplitude? Can the dynamics of the spinal circuits be adjusted using non-specific central control signals according to a simple scaling rule? Or does the CNS have to learn an explicit mapping between desired kinematics and control parameters?

We have speculated in the introduction that it might be intersegmental force feedback carried by Ib afferents that provides the mechanism by which the nervous system compensates for interaction torques. Such signals reliably encode muscle force (Mileusnic and Loeb, 2009), can be adjusted in sensitivity through interaction with Ia afferents (McCrea, 1992), and result in widespread modulation of motor neurons innervating muscles acting at adjacent joints (Jankowska et al., 1981). Moreover, functionally deafferented patients in some cases make systematic movement errors indicative of a failure to counteract interaction forces, which demonstrates a functional role for proprioception in the compensation of internal loads (Ghez and Sainburg, 1995; Sainburg et al., 1995). And motion-dependent feedback across spinal segments has been shown to modulate ongoing limb dynamics in the cat (Smith and Zernicke, 1987; Koshland and Smith, 1989). Our model only allows us to confirm that force feedback can indeed play a role in the compensation of internal loads, though we have not identified precisely what that role is. It is clear from our results that without the contribution of Ib afferents the performance of the model is greatly diminished. But both intra- as well as intersegmental effects of Ib afferents seem to contribute to the appropriate modulation of spinal neurodynamics. Further work is required to separate and study these two effects in more detail.

Our model does not address the question of how spinal circuitry might come to be organized in the manner presented here. The evolutionary optimization procedure does not serve as a model of how the appropriate connectivity could be learned. Nevertheless, it is known that neural circuits are subject to activity-dependent plasticity both in the developing as well as the mature spinal cord (Changeux and Danchin, 1976; Nelson et al., 1990; Lo and Poo, 1991; Wolpaw and Carp, 1993; Schouenborg, 2003; Wolpaw, 2010; Tahayori and Koceja, 2012). Moreover, the reciprocal structure typical of spinal circuits can arise through self-organization enabled by simple Hebbian-like learning rules in initially undifferentiated networks undergoing spontaneous activity (van Heijst et al., 1998; Petersson et al., 2003; Marques et al., 2013). Together with models demonstrating the feasibility of trial-and-error learning (Raphael et al., 2010) this evidence suggests that acquisition of appropriately tuned neural circuits in the spinal cord is possible.

Further work is needed to investigate the exact role of the different feedback modalities (position, velocity, and force) and interneurons in the accommodation of interaction torques in model spinal circuits. Given that the organization of the human spinal cord in reality is significantly more complicated than modeled here, and our substantial yet incomplete knowledge regarding its structure and functionality, an interesting avenue to explore would be to study the complete ensemble of spinal models that fill in unknown data and conform with behavioral and neurophysiological data. Such a methodology has been applied, for example, to study the mechanism underlying klinotaxis in *C. elegans* and to propose experiments that can distinguish between different hypotheses regarding its neural implementation (Izquierdo and Beer, 2013).

In conclusion, the work presented here demonstrates the feasibility of equilibrium point control for multijoint reaching movements subject to varying intersegmental loads. The model shows that internal models and predictive compensation of such loads are not required for the range of movements studied here. But it does not allow us to refute the possibility that such mechanisms are indeed used by the CNS for this or other purposes. The model also indicates that EP style control of reaching movements is dependent on a mapping of desired movement kinematics to control parameters and on the appropriate self-organization of spinal circuitry. We do not propose that it is the sole, or even most central, function of spinal circuitry to implement equilibrium point control, or to coordinate different muscles for interaction torque accommodation. Indeed, while other models not incorporating spinal interneurons might also be able to accommodate interaction torques to some extent, reflex circuitry in the spinal cord in conjunction with central modulation may allow for greater flexibility in the execution of a given class of movements (such as adaptation to varying energy, speed and accuracy trade-offs), and may underlie the ability to perform different classes of mechanical action (such as control of position, force or stiffness) as and when needed.

# **FUNDING**

This work is funded by the project "eSMCs: Extending Sensorimotor Contingencies to Cognition," FP7-ICT-2009-6 no: 270212.

#### **ACKNOWLEDGMENT**

The authors are grateful to Marieke Rohde, Hugo Gravato Marques and Mike Beaton, as well as the reviewers of this manuscript, for their comments and recommendations on an earlier version of this article.

# **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fncom. 2014.00144/abstract

## **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 30 May 2014; accepted: 23 October 2014; published online: 11 November 2014.*

*Citation: Buhrmann T and Di Paolo EA (2014) Spinal circuits can accommodate interaction torques during multijoint limb movements. Front. Comput. Neurosci. 8:144. doi: 10.3389/fncom.2014.00144*

*This article was submitted to the journal Frontiers in Computational Neuroscience.*

*Copyright © 2014 Buhrmann and Di Paolo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Emulated muscle spindle and spiking afferents validates VLSI neuromorphic hardware as a testbed for sensorimotor function and disease

# *Chuanxin M. Niu1, Sirish K. Nandyala2 and Terence D. Sanger 2,3,4\**

*1 Department of Rehabilitation, Ruijin Hospital, School of Medicine, Shanghai Jiao Tong University, Shanghai, China*

*2 Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, USA*

*3 Biokinesiology, University of Southern California, Los Angeles, CA, USA*

*4 Neurology, University of Southern California, Los Angeles, CA, USA*

#### *Edited by:*

*Misha Tsodyks, Weizmann Institute of Science, Israel*

*Reviewed by: Eric Jon Perreault, Northwestern University, USA Richard Nichols, Georgia Tech, USA*

#### *\*Correspondence:*

*Terence D. Sanger, Sanger Lab, Department of Biomedical Engineering, University of Southern California, 1042 Downey Way, Los Angeles, CA 90089, USA e-mail<EMAIL_ADDRESS>
The lack of multi-scale empirical measurements (e.g., recording simultaneously from neurons, muscles, whole body, etc.) complicates understanding of sensorimotor function in humans. This is particularly true for the understanding of development during childhood, which requires evaluation of measurements over many years. We have developed a synthetic platform for emulating multi-scale activity of the vertebrate sensorimotor system. Our design benefits from Very Large Scale Integrated-circuit (VLSI) technology to provide considerable scalability and high-speed, as much as 365× faster than real-time. An essential component of our design is the proprioceptive sensor, or muscle spindle. Here we demonstrate an accurate and extremely fast emulation of a muscle spindle and its spiking afferents, which are computationally expensive but fundamental for reflex functions. We implemented a well-known rate-based model of the spindle (Mileusnic et al., 2006) and a simplified spiking sensory neuron model using the Izhikevich approximation to the Hodgkin–Huxley model. The resulting behavior of our afferent sensory system is qualitatively compatible with classic cat soleus recording (Crowe and Matthews, 1964b; Matthews, 1964, 1972). Our results suggest that this simplified structure of the spindle and afferent neuron is sufficient to produce physiologically-realistic behavior. The VLSI technology allows us to accelerate this behavior beyond 365× real-time. Our goal is to use this testbed for predicting years of disease progression with only a few days of emulation. This is the first hardware emulation of the spindle afferent system, and it may have application not only for emulation of human health and disease, but also for the construction of compliant neuromorphic robotic systems.

**Keywords: spindle, motor control, afferents, emulation, neuromorphic**

# **INTRODUCTION**

The multi-scale nature of human nervous system makes it difficult to measure all relevant information about sensorimotor function. Moreover, such inadequacy of multi-scale information also prevents us from understanding the mechanism of neurological diseases, whose expression depends on networks and interconnectivity of neurons. In the case of movement disorders in childhood, the ultimate clinical effect of cellular injury may take years or even decades to fully emerge (Sanger, 2003) due to a complex interplay among the child's injury, development, and experience. Therefore, the lack of detailed and complete sensorimotor measurements, both in multi-scale and across long timespan, makes it extremely difficult to quantify the mechanism of disease or to predict the efficacy of long-term treatments.

Here, we attempt to understand sensorimotor functions using a synthetic approach. This requires creating sensorimotor models that are biologically realistic but emulated on nonbiological hardware. The history of synthetic approaches in neuroscience dates back to the 1940s, when scientists started creating artificial neurons and neural networks using electronic circuits (McCulloch and Pitts, 1943). Models of neuron dynamics (Hodgkin and Huxley, 1952; Rosenblatt, 1958) soon emerged to be simulated on digital computers. Since the 1980s, specialpurpose hardware with Very Large Scale Integrated-circuit (VLSI) technology started to benefit from some of the key insights in neural computation, including asynchrony among neurons, spike representation of information, and self-improving mechanisms such as plasticity (Mead, 1989; Serrano-Gotarredona et al., 2009; Indiveri et al., 2011; Neftci et al., 2014). This category of designs, termed "neuromorphic" hardware, has been successful in understanding mechanisms of memory (Chicca et al., 2003), visual representation (Lichtsteiner et al., 2008), and recently cognitive function (Eliasmith et al., 2012). For sensorimotor function, the synthetic approach should not only describe neurons, but also the physiological environment (muscles, proprioceptors, peripheral nerves, skeletal system) that the neurons interact with. If the goal is to predict long-term changes of sensorimotor function, the emulation speed must exceed real-time. These requirements pose a major challenge when designing the hardware, as the running speed of may decrease significantly with model complexity. In this case, our approach for sensorimotor modeling must accommodate enough detail for multi-scale simulation, but it must also maintain high speed when scaling to larger network sizes.

Field Programmable Gate Array (FPGA), a programmable VLSI device, can parallelize neuron simulations at very high speed (Guerrero-Rivera et al., 2006; Li et al., 2010). The inherent parallel structure of FPGAs permits that when model complexity increases, the running speed can be maintained by adding units. The parallel design uses each clock cycle to update multiple neurons, access a wide range of memory, and send out neuron spikes for inter-chip communication; therefore FPGAs are able to accomplish more calculations than clock-cycled computers at the same central clock frequency. In our previous work (Niu et al., 2012), we showed the technical details of the emulation platform using FPGAs, which achieved a running speed 365× faster than real-time for a simulation of the spinal components of the monosynaptic stretch reflex. In the current study, we simulate the muscle spindle proprioceptors and their attached sensory afferent neurons. The spindle system is an essential part of the monosynaptic stretch reflex, as well as a major source of state information for feedback control and stabilization. Because of the model complexity of spindles, there has not previously been an attempt to simulate them in hardware. Therefore, our high-speed hardware simulation of the muscle spindle will provide an essential element for studying the stretch reflex as well as for understanding stabilization and feedback control of biological and neuromorphic motor systems.

Among the many studies investigating muscle spindles in mammals (Eldred et al., 1962; Lennerstrand, 1968; Boyd et al., 1977; Hulliger, 1984), the series of experiments conducted by Matthews and colleagues (Crowe and Matthews, 1964b; Matthews, 1964, 1972) characterized cat soleus spindles in considerable detail. Their representative data provide the reference for the rate-based spindle model (Mileusnic et al., 2006) chosen for our study. We expect our spike-based emulation to enrich the original rate-based spindle model: all spiking afferents will be available for analysis, and the spiking behavior should be compatible with physiological recording. The advantage of converting to a spike-based representation is that we can compare our results directly against physiologically-measured signals, we can investigate questions related to the number or firing pattern of sensory afferents, and we can investigate the effect of spike timing on plasticity. For the purpose of emulating sensorimotor function and disease, we argue that the model of spindle afferents does not have to include all anatomical details, but its outcome must satisfy at least two constraints: (1) distinguishable firing pattern between Group Ia and Group II afferents, and the difference should be compatible with physiological data; (2) distinguishable change in behavior that reflects changes in dynamic and static gamma fusimotor drive. We verified the VLSI emulation against these two constraints. Because of the physiological differences between cats and humans (Prochazka and Hulliger, 1983), we also directly compared our emulated signals with human spindle afferents (Edin and Vallbo, 1990) to see whether this model can be used to emulate human data. If successful, our results will provide a testbed with sufficient detail to represent known sensorimotor physiology in order to describe and predict the effect of neural injuries over many years.

# **MATERIALS AND METHODS**

#### **APPARATUS**

The emulated neurological system comprises a muscle spindle with three types of fibers (bag1, bag2, and chain), two groups of afferents (Group Ia and II) responding to two types of fusimotor drives (gamma dynamic and static), illustrated in **Figure 1A**. The mathematical descriptions of muscle spindle are implemented using (FPGA, Xilinx Spartan-6), a programmable version of VLSI electronic chips.

We favor FPGAs over pipelined hardware such as Graphic Processing Units (GPUs) or clustered CPUs due to its inherent parallelism that resembles neural circuitry. It is also argued that when networking multiple units for large-scale disease emulation, FPGAs are significantly more flexible for enabling custom-built communication protocols, including neuromorphic transmission protocols that directly transmit neuron-like spikes.

#### **FPGA HARDWARE DESIGN**

The implementation requires translating the equations of neurons and spindles into digital circuits described in Verilog Hardware Description Language (HDL). We enforce a combinatorial circuit design as opposed to a sequential clocked design when translating all models, therefore we can maximize performance since circuits do not have to wait for clock signal boundaries. A central scheduler periodically updates all emulated models and distributes the latest information to corresponding parts. Each update generates neurological activities accounting for 1 ms in the real world. The frequency of update can be easily adjusted on digital VLSI, therefore when increasing the updating frequency the emulation will operate faster than real-time. The speed of emulation can keep increasing as long as all models finish their state update within each accelerated clock cycle. The maximal speed of emulation is constrained by how fast the electronic signals can propagate through the combinatorial circuits.

**Figure 1B** shows a working setup of the hardware platform. The spiking activities of up to 80 neurons are measurable directly from the FPGA using an oscilloscope for "virtual neuron recording," which on-site verifies the accuracy of emulation with hardware in the loop. Other activities are transferred to a data-logging computer for offline analysis. Note that when collecting data for offline analysis, the emulation has to be slowed below real-time due to limited bandwidth between the FPGA and the datalogging computer; the emulation is accelerated to 365× real-time only when investigating long-term changes. The FPGA communicates with the data-logging computer through OpalKelly development kits (XEM6010, OpalKelly Inc.). More technical details can be found in Niu et al. (2012).

#### *Floating-point arithmetics in combinational logic*

Spindle models are evaluated in IEEE-754 single-precision floating-point numbers. Typical floating-point arithmetic IP

cores are either pipe-lined or based on iterative algorithms such as CORDIC, all of which require clocks to schedule the calculation. In our testbed, no clock is allowed for model evaluation thus all arithmetics need to be executed in pure combinational logic. The implementations of adder and multiplier are inspired by the open source project "Free Floating-Point Madness," available at http://www.pldworld.com/\_hdl/2/\_ip/www.hmc.edu/ chips/fpmul.html. The modified code used in this study is available upon request.

Floating-point division is more resource demanding than multiplications. We approximated the division with additions and multiplications therefore avoided the direct implementation of floating-point division. Our approach is inspired by an algorithm described by Lomont (2003), which provides a good approximation of the inverse square root for any positive number *x* within one Newton-Raphson iteration:

$$Q(\mathbf{x}) = \frac{1}{\sqrt{\mathbf{x}}} \approx \mathbf{x} \left( 1.5 - \frac{\mathbf{x}}{2} \cdot \mathbf{x}^2 \right) (\mathbf{x} > \mathbf{0}) \tag{1}$$

where Q (x) only contains additions and multiplications. Any division with a positive divisor can be achieved if two blocks of Q (x) are concatenated:

$$\frac{a}{b} = \frac{a}{\sqrt{b} \cdot \sqrt{b}} = a \cdot Q(b) \cdot Q(b) \,\mathrm{(b > 0)}\tag{2}$$

It is trivial to adjust this algorithm for negative divisors (b < 0).

#### *Serialize neuron evaluation using time-shared multiplexing*

Consider that the spindle model only needs to update at 365 kHz even in the highest speed (1 ms time granularity, 365× real-time), there is room for time-sharing the FPGA logic-gates among neurons, thus using fewer logic-gates to emulate larger amount of neurons. The maximal number of neurons that can be serialized (*Nserial*) is constrained by the following relationship:

$$C \times N\_{serial} \times \text{365} \times F\_{emu} \le F\_{f\text{çça}} \tag{3}$$

Here *Ffpga* is the fastest clock rate that a FPGA can operate on; *C* is the minimal clock cycles needed for updating each state variable in the on-chip memory, in our case *C* = 2 due to an optimized design for memory access; *Femu* = 1 *kHz* is the time granularity of emulation (1 ms), and 365 × *Femu* represents 365× real-time. Consider that Xilinx Spartan-6 FPGA devices have a maximum 200 MHz central clock frequency, the theoretical maximum of neurons that can be serialized is

$$N\_{serial} \le \frac{200 \, MHz}{2 \times 365 \times 1 \, kHz} \approx 274 \tag{4}$$

In the current design we choose *Nserial* = 128.

#### **EMULATION OF MUSCLE SPINDLES**

As the sensory organ that provides the main source of proprioceptive information, a typical muscle spindle produces afferents Group Ia (transducing muscle length and a function of velocity) and Group II (primarily sensing muscle length). The activity of spindle is modulated by two types of fusimotor drive (gamma dynamic and gamma static). At present there are no analytical equations that can accurately describe the dynamics of a spindle. However, spindle behavior can be approximated using ordinary differential equations for numerical solutions. Though few such models are available, the one presented by Mileusnic et al. (2006) showed a close fit to the firing rate recorded from cat soleus spindles. We implemented this spindle model as the first step to deduce the essential elements for realistic spindle afferents.

The chosen spindle model used differential equations to describe the relationship between afferent firing rates and the instantaneous muscle status (length and lengthening velocity). The overall structure of the spindle model can be summarized as follows:

$$\dot{R}(t) = f\left(L, \dot{L}, \,\chi\_{\text{dynamic}}, \,\varphi\_{\text{static}}, t\right) \tag{5}$$

$$
\hat{R}\left(t\right) = \int \dot{\vec{R}}\left(t\right)dt\tag{6}
$$

where *R* (*t*) is the instantaneous firing rate of an afferent fiber at time *t*, *L* is the length of the muscle, γ*dynamic* and γ*static* are the firing rates of gamma fusimotor drives. This model described each spindle fiber as coupled springs and damping; it characterized how three types of fiber (bag1, bag2, nuclear chain) contribute to each group of spindle afferent (Ia and II). The thixotropic property of spindle, i.e., that the spindle output depends on the history of stretching (Hasan and Houk, 1975), was not captured in this model.

A full combinatorial implementation of the spindle model requires more silicon resources than available on a single Spartan-6 FPGA chip. Therefore, we optimized the model by trimming its transcendental functions and replacing some on-chip calculations with pre-calculated values. For example, a subset of the model for bag1 fiber is shown below:

$$\dot{\mathbf{x}}\_0 = \left(\frac{\Gamma\_{dynamic}}{\Gamma\_{dynamic} + \Omega^2} - \mathbf{x}\_0\right) / \mathbf{r} \tag{7}$$

$$
\dot{\mathbf{x}}\_1 = \mathbf{x}\_2 \tag{8}
$$

$$\dot{\mathbf{x}}\_2 = \frac{1}{M} \left[ T\_{SR} - T\_B - T\_{PT} - \Gamma\_1 \mathbf{x}\_0 \right] \tag{9}$$

where

$$T\_{\rm SR} = K\_{\rm SR} \left( L - \mathbf{x}\_{\rm I} - L\_{\rm SR0} \right) \tag{10}$$

$$T\_{\rm SR} = (B\_0 + B\_1 \chi\_0) \cdot (\chi\_1 - R) \cdot \text{CSS} \cdot |\chi\_2|^{0.3} \tag{11}$$

$$T\_{PR} = K\_P R \left(\mathbf{x}\_1 - L\_P R \mathbf{0}\right) \tag{12}$$

refer to Mileusnic et al. (2006) for the definitions of variables. As can be seen from Equation (11), the absolute value of lengthening velocity (*x*2) was raised to the power of 0.3, which is trivial to program in Matlab/C++ but has no corresponding syntax in Verilog. We approximated |*x*2| 0.3 as |*x*2| 0.25, because the power of 0.25 is straightforward to acquire by iterating on *Q* (*x*) of Equation (1):

$$Q\left(Q\left(\mathbf{x}\right)\right) = \left(\left(\mathbf{x}^{-0.5}\right)\right)^{-0.5} = \mathbf{x}^{0.25}\left(\mathbf{x} > 0\right) \tag{13}$$

In addition, most of the equations were expanded to polynomials, such that the constant parts in the equations were pre-calculated without consuming multipliers, e.g., 1/*M* in Equation (9) was pre-calculated as a single number instead of doing an on-chip division.

After optimization we are able to host an entire spindle using only 60% of a single Spartan-6 FPGA chip, with a running speed peaking at 500× real-time. Typically a single Spartan-6 chip can support one spindle connected with at least 1024 neurons at 365× faster than real-time.

#### **EMULATION OF SPIKING NEURONS**

From the perspective of information theory, the function of a neuron is to convert post-synaptic currents into a train of binary spikes with limited bandwidth (Sanger, 2011). In an emulated neurological system focusing on its functional role, a neuron can be modeled to any level of detail as long as it satisfies the protocol of "current-in, spike-out." In the current study we adopted the neuron model developed by Izhikevich (2003), which approximates the Hodgkin–Huxley model (Hodgkin and Huxley, 1952). In our emulation, we set the four parameters (a, b, c, d) required in the Izhikevich model to ensure all sensory neurons fit Hodgkin's description of Class 1 excitatory neurons (Hodgkin, 1948; Izhikevich, 2007). Class 1 excitatory neurons are one of the major types found in human spinal cord recordings (Prescott et al., 2008). Since the firing rate of Class 1 excitatory neurons is a monotonic representation of post-synaptic current over a large dynamic range, it allows straightforward conversion from the spindle output (in firing rate) to the input of the neuron (in post-synaptic current).

Although we enforce a pure combinatorial design to maximize the running speed, neurons still time-share physical circuits in order to maximize the population size with limited silicon surface. Pseudorandom white uniform noise (5 mV amplitude) was added to the membrane potential of each neuron to introduce variability for population firing. The noise is introduced to represent the large number of inputs that a neuron usually receives from the dendritic tree. The noise level was set to create a typical 4.8 mV fluctuation in the membrane potential (Fellous et al., 2003). Pseudorandom noise is generated using a linear-feedback-shift-register (George and Alfke, 2001).

#### **SPIKING RESPONSES OF SPINDLE TO VIRTUAL MUSCLE STRETCH**

There has been a long-lasting history studying rapid excitatory responses of a muscle following stretch that dates back to 1751 by Robert Whytt (Pearce, 1997). Over the past centuries it has been revealed that the process of response is a complex muscle reaction, with multiple excitatory responses occurring at different latencies following a muscle stretch (Hammond, 1955). We emulated a series of classic muscle stretch experiments that were originally performed in cat soleus muscles (Crowe and Matthews, 1964b; Matthews, 1964, 1972). Since these datasets are matched in the original spindle model (Mileusnic et al., 2006) to produce compatible firing rates, our emulation is intended to produce compatible spike patterns under similar muscle stretch. Four types of stretch stimulus were introduced to the emulated spindle: linear stretch, tap, sinusoidal stretch and release. We focus on whether the spindle afferents show: (1) distinguishable firing pattern between Group Ia and Group II afferents, and the difference should be compatible with physiological data; (2) distinguishable change in behavior that reflects changes in dynamic and static gamma fusimotor drive.

We further analyzed the emulated spindle activity using a white-noise approach, which has been widely used in system identification for non-linear biological systems (Marmarelis, 2004). We stimulated the emulated spindle using low-pass filtered white noise. It is expected that the firing rate reconstructed from our emulation should be statistically correlated with the firing rate produced by the rate-based spindle model. In addition, the emulated spindle activity was compared to spindle afferent recordings from humans. Due to the known difference between spindles of cats and humans (Prochazka and Hulliger, 1983), we do not expect our emulation to exactly match human data but the comparison should reflect the known difference.

It is noteworthy that all virtual recordings in this study required slow emulation below real-time, otherwise the amount of data produced by FPGA per unit time would exceed the bandwidth for data logging. The mathematical correctness of emulation when accelerated to the full 365× real-time speed is ensured by in-loop fault checks.

## **RESULTS**

#### **QUALITATIVE DIFFERENCE BETWEEN Ia AND II AFFERENT SPIKES**

We first replicated the stretch-and-hold experiment (Matthews, 1972) conducted with the soleus muscle of a decerebrate cat. In their work, the cat soleus muscle was stretched by 14 mm within 200 ms with the muscle maintained elongated after stretching. The spiking responses of cat were measured using intramuscular recording and are shown in **Figure 2A** (reproduced with permission). Responses were recorded both in the presence (ventral roots intact) and absence (ventral roots cut) of tonic fusimotor activity. In both cases the Primary (Ia) afferent exhibited stronger phasic response than Secondary (II) afferent, especially when the stretching velocity is greater than zero. In contrast, the Secondary (II) afferent produced stronger tonic response during the entire process of stretching.

Equivalent stretch-and-hold experiments were tested in VLSI emulation. Note that the original spindle model requires all muscle lengths normalized to the length of relaxed muscle (resting length, *L*0). The spindle was stretched to 36.8% of *L*0, corresponding to the 14 mm elongation of cat soleus muscle (Matthews, 1972) from an average rest length of 38 mm (Scott et al., 1996). **Figure 2B** shows the emulated action potentials in response to the 0.368*L*0 virtual stretch lasting for 200 ms. We observed a clear distinction between Ia and II afferents in the emulated results. It can also be seen that in the presence of gamma fusimotor drive (ventral roots intact), both Ia and II afferents were more active in a similar way between real and emulated responses. The difference between Ia and II afferents is qualitatively similar to physiological data.

Note that in **Figure 2A** the muscle tensions were presented to show the muscle's mechanical response to external stretch, whereas the actual waveform of stretch should be in the shape shown in **Figure 2B** (muscle length). In our emulation, the muscle tension was only implicitly calculated for the spindle as an intermediate variable, and was hence not directly measurable.

#### **QUALITATIVE DIFFERENCE BETWEEN STATIC AND DYNAMIC FUSIMOTOR DRIVE**

We also compared the biological and emulated afferents from a primary afferent neuron, when the spindle fiber is stimulated with different fusimotor drives. In particular, two fusimotor drives (gamma dynamic and static) were selectively stimulated during a sinusoidal stretch of 1 mm peak-to-peak at 3 Hz (Crowe and Matthews, 1964a). The original recording is shown in **Figure 3A**. External stimuli to the gamma fibers were activated between the 3rd and 4th cycle (horizontal line in **Figure 3A**). As expected, the effect of gamma dynamic drive was to facilitate the phasic response of cat soleus spindle, while the tonic responses were mostly facilitated by gamma static drive. Emulated spindle afferent spiking is shown in **Figure 3B**. A virtual sinusoidal stretch of 0.026*L*0 peak-to-peak at 3 Hz was applied to the emulated spindle. Gamma fusimotor drives were initially set to 0 Hz and activated to 80 Hz between the 3rd and 4th cycle. Again we observed distinctive spiking patterns in response to gamma static and dynamic stimuli.

#### **POPULATION FIRING IN RESPONSE TO VARIOUS MUSCLE STRETCHES**

One significant advantage of neuromorphic emulation is that the spikes of a large number of neurons can be stored and analyzed together with many other emulated activities. Here we test whether the neuron ensemble shows expected spiking behaviors in response to various types of stretching waveforms. **Figure 4A** shows four types of stretch waveforms (linear stretch, tap, sinusoidal stretch and release) along with hand-drawn schematized spike responses summarized from experimental observations (Matthews, 1964) in absence of fusimotor drive, reprinted with permission.

The emulated responses are shown in **Figure 4B**. In all cases except for muscle release, the emulated spike train from a sample neuron matches the empirical responses. In addition, both the Primary and Secondary ensembles, each comprising 128 neurons, show variable but congruent patterns compared to a single neuron. The only noticeable exception is in muscle release (**Figure 4B**, rightmost column), where the Primary afferents show a burst of firing after the release stops. This burst is caused by the fact that a sudden stop is equivalent to a momentary stretch at very high velocity. It was not schematized in **Figure 4A**, however, probably because the original figure was hand-drawn such that some details were unnecessarily omitted. According to Poliakov and Miles (1994), sudden bursts of electromyography (EMG) were indeed present at the end of masseter releasing. Such discrepancy demonstrated the ability of our platform to validate inconsistent observations.

representative neuron and the neuron ensemble when the four types of muscle stretching were applied. The raster shows 128 Primary afferent neurons and 128 Secondary afferent neurons. Changes in spike pattern are similar between experimental and emulated results, with emulated results showing irregularities in spiking due to added random noise. Note that in the case of muscle release (last column), the Primary afferents show a burst of firing after the release suddenly stopped. This burst was not schematized in **(A)**, probably because sudden stops are usually smoothed significantly by skin, ligaments, and other tissues.

In addition, there is a noticeable distinction between the regularity of spikes of **Figures 4A,B**. One explanation is that our emulation explicitly introduced synaptic noise, and therefore the individual neuron will exhibit a Poisson-like random firing pattern. This particular randomness was perhaps not the focus when **Figure 4A** was schematized.

Another feature of muscle spindle is that its afferents usually spike at a significantly high rate during the beginning of muscle stretch, resulting in an "initial burst." This feature was measure in rats by Haftel et al. (2004) as shown in **Figure 5A**. We used the same triangular waveform (**Figure 5B**) to stretch our emulated spindle for 25 repetitions. The instantaneous frequency averaged across 25 repetitions (**Figure 5C**) showed similar initial bursts compared to experimental data.

#### **CORRELATION BETWEEN RATE-BASED MODEL AND SPIKE-BASED EMULATION**

We tested the spindle afferents using a pseudo-white-noise waveform, which examines the input-output relationship of a dynamic system with rich frequency components (Marmarelis, 2004). The waveform includes a series of pseudorandom numbers low-pass filtered at 5 Hz cutoff frequency. The low-pass filtering screens out the unrealistic high-frequency stretching that is usually damped by skin and ligaments. One sample of a 4-s stretch is shown

in **Figure 6A**. The corresponding Primary Ia response in firing rate (**Figure 6B**) generates stochastic spiking shown in a raster plot (**Figure 6C**, raster of 16 neurons). The firing rate produced by the original spindle model is compared to the total spike count recorded from the 256 emulated spindle afferent neurons (**Figure 6D**). The total spike count shows irregularity and variability in each run because of the stochastic spiking in the emulated neuron population.

We stretched the spindle for an equivalent of 160 s in realtime using the low-pass filtered pseudo-white-noise waveform. The total spike count is significantly correlated with the expected firing rate for both Primary Ia (*p* < 0.0001, *r* = 0.813) and Secondary II (*p* < 0.0001, *r* = 0.810). These results verify that our emulation provides a consistent spike-based representation of muscle length and lengthening velocity; the spiking outcome is statistically compatible with the original rate-based spindle model on which it is based.

#### **COMPARISON TO HUMAN SPINDLE AFFERENTS**

Although the original spindle model was developed based on cat soleus spindles, we replicated the experiments done with human spindle to compare our emulation with human data. We introduced 0.7*L*0 stretch using the waveform reported by Edin and Vallbo (1990). The human spindle recording and emulated activity are shown in **Figure 7**.

As can be seen, the overall spike rate is lower in emulation than human data; the emulated spike rates show less contrast between long and short muscle length as compared to human data, i.e., reduced static response. These differences are compatible with the finding that human spindles contain more intrafusal fibers, which are probably dominated by nuclear chain (static) fiber (reviewed

in Prochazka and Hulliger, 1983). Quantitative comparison is not performed due to limited human data. In principle, the spindle model could be tuned to fit human spindles by manipulating the damping terms of the spindle fiber.

#### **DISCUSSION**

Using our recently developed technique of neuromorphic emulation on programmable digital VLSI hardware, we showed that a collection of spiking afferent fibers driven by a detailed model of muscle spindles suffice to produce biologically-realistic spindle afferents. The observed firing pattern from emulated spindle neurons matches classical intramuscular recordings from cat soleus muscle, and the emulated responses to gamma fusimotor drive show no qualitative difference from experimental recordings. All emulations can be accomplished at 365× real-time, which allows estimating long-term changes and large population behavior efficiently. The emulated spindles differ from human spindle activity but the difference is compatible with the known difference between cat and human spindles. Our results provide a strong validation of using neuromorphic emulation as a testbed for neurophysiological studies. We can now test the roles of more complex structures, such as realistic muscles, inhibitory neurons, or supra-spinal circuitry in producing movement behavior. The multi-scale design enables us to emulate pathological conditions that are physiologically tenable, e.g., death of neurons or absence of gamma fusimotor drive. It creates an essential step toward investigating how such pathological conditions could contribute to disease progression in childhood.

Synthetic emulation using neuromorphic hardware provides three major advantages compared to empirically studying the real biological system. First, it isolates the subsystem-of-interest from the compounding factors that are very difficult to tease out *in-vivo*, including studying the spinal reflex isolated from

supra-spinal influence. Second, it allows interventions that are usually difficult to introduce when studying biological systems, e.g., differentially adjusting the relative weight of gamma dynamic vs. static drive. Third, the hardware acceleration allows for predicting the system's long-term development with significantly less time. Our vision for the synthetic approach is not to replace empirical studies but rather to inspire testable hypotheses for experiments and clinical applications. Another important feature is that this platform can be used to verify the source of motor variability from different physiological origins, such as the intrinsic firing variability of neurons, motor noise associated with muscles, or the individual properties of mechanoreceptors.

The purpose of this study is to validate the neuromorphic hardware as a testbed for spindle activity, therefore we focused on (1) implementing a selected model of spindle (Mileusnic et al., 2006), (2) adapting it to enable spiking afferents, (3) accelerating it to 365× real-time; we did not focus on improving the model to match more experimental data than it originally could. Nevertheless, the limitations of this spindle model should be acknowledged. One major limitation is that the longknown thixotropic property of spindle activity (Hasan and Houk, 1975) was not captured; also the model was based on cat data, thus it must be re-calibrated for modeling human movement disorders. It is worth noting that other excellent models of spindle also exist (e.g., Hasan, 1983; Lin and Crago, 2002), they either used a less computationally expensive approach to model the non-linear velocity dependency (Hasan, 1983), or succeeded in unifying spindle and Golgi Tendon Organ with the same structure (Lin and Crago, 2002). Recent work also refined the original spindle model by examining the non-linearity in its components (Lan and He, 2012). Our design of the platform is open and flexible for including additional features in future improvements, or switching to different models if necessary.

The original model of spindle (Mileusnic et al., 2006) focused on the firing rates of spindle afferents instead of their spiking patterns. Such rate-based models are incompatible with our overall goal of disease emulation, where the spike timing is crucial for motor variability and neuroplasticity. In acknowledgement of the soundness of the original rate-based spindle model, our major improvement is to enable a large number of spiking neurons driven by a spindle. Moreover, the original model was developed in Matlab Simulink operating slower than real-time, while our hardware implementation permits faster-than-real-time performance. This is the first physical, portable, and realistic proprioceptor that can provide synthetic proprioception for robots, virtual neurophysiological studies, and prediction of clinical outcomes.

#### **ACKNOWLEDGMENTS**

The authors are grateful for support from the National Institute of Neurologic Disorders and Stroke (R01-NS069214), and the James S. McDonnell Foundation. The authors thank Dr. Gerald Loeb for the help with the spindle model.

#### **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 08 June 2014; accepted: 17 October 2014; published online: 04 December 2014.*

*Citation: Niu CM, Nandyala SK and Sanger TD (2014) Emulated muscle spindle and spiking afferents validates VLSI neuromorphic hardware as a testbed for sensorimotor function and disease. Front. Comput. Neurosci. 8:141. doi: 10.3389/fncom.2014.00141 This article was submitted to the journal Frontiers in Computational Neuroscience.*

*Copyright © 2014 Niu, Nandyala and Sanger. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Hybrid model of the context dependent vestibulo-ocular reflex: implications for vergence-version interactions

# *Mina Ranjbaran\* and Henrietta L. Galiana*

*Department of Biomedical Engineering, McGill University, Montreal, QC, Canada*

#### *Edited by:*

*Misha Tsodyks, Weizmann Institute of Science, Israel*

#### *Reviewed by:*

*Petia D. Koprinkova-Hristova, The Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, Bulgaria Arthur Prochazka, University of Alberta, Canada*

#### *\*Correspondence:*

*Mina Ranjbaran, Department of Biomedical Engineering, McGill University, 3775 University, Montreal, QC H3A2B4, Canada e-mail: mina.ranjbaranhesar<EMAIL_ADDRESS>
The vestibulo-ocular reflex (VOR) is an involuntary eye movement evoked by head movements. It is also influenced by viewing distance. This paper presents a hybrid nonlinear bilateral model for the horizontal angular vestibulo-ocular reflex (AVOR) in the dark. The model is based on known interconnections between saccadic burst circuits in the brainstem and ocular premotor areas in the vestibular nuclei during fast and slow phase intervals of nystagmus. We implemented a viable switching strategy for the timing of nystagmus events to allow emulation of real nystagmus data. The performance of the hybrid model is evaluated with simulations, and results are consistent with experimental observations. The hybrid model replicates realistic AVOR nystagmus patterns during sinusoidal or step head rotations in the dark and during interactions with vergence, e.g., fixation distance. By simply assigning proper nonlinear neural computations at the premotor level, the model replicates all reported experimental observations. This work sheds light on potential underlying neural mechanisms driving the context dependent AVOR and explains contradictory results in the literature. Moreover, context-dependent behaviors in more complex motor systems could also rely on local nonlinear neural computations.

**Keywords: sensory-motor mapping, vestibulo-ocular reflex, context dependent reflex, mathematical model, disconjugate eye movement, ocular nystagmus**

# **1. INTRODUCTION**

The vestibulo-ocular reflex is an involuntary eye movement that stabilizes gaze in space during head movements for clear and blur-free vision. The rather simple neural substrate of the VOR, the so-called *three neuron arc* (de No, 1933), makes it an appropriate model to study sensory-motor behavior. Rotational and translational head movements are sensed by the vestibular system (semicircular canals and the otolith organs) in the inner ear. Vestibular afferents relay sensory information to the vestibular nuclei (VN) and prepositus hypoglossi (PH) centers in the brainstem. These centers act as the main controller and combine sensory signals with internal efference copies of the controlled plant(s), eye orientation, to drive motor-neurons appropriately. Extraocular muscles then apply torques on the eyeball that result in the eye movements.

VOR nystagmus consists of compensatory (slow phase) and reorienting (fast phase) segments. The slow phases of the VOR stabilize gaze in space by moving the eyes in the opposite direction to the head movement, while the fast phases redirect the gaze at high speeds in the direction of the head velocity. We focus on the angular VOR (AVOR), tested with passive whole-body rotation in the dark while recording conjugate (eye movements in the same direction) or monocular horizontal eye movements. **Figure 1** shows an example of the VOR during sinusoidal head rotations using electrooculography (EOG) in the dark: some slow and fast phase segments are marked. The sawtooth-like pattern of the eye movement is a characteristic of most types of eye movements and is known as ocular nystagmus. In clinical tests, the VOR is characterized by its *gain* defined as the ratio of peak eye velocity to peak head velocity during harmonic testing or short pulse perturbations.

While the head movements initiate the VOR, this reflex is also influenced by contextual factors such as viewing distance (Viirre et al., 1986; Crane and Demer, 1998). Since the eyes are not centered on the head, holding gaze on a near target requires more ocular rotation than for a relatively far target during head movements. In other words, the AVOR gain increases as a function of decreasing fixation distance, that can be described geometrically. The majority of models that attempted to explain target-distance dependent VOR responses relate this property to (i) an internal signal proportional to the inverse of target distance that scales VOR gain (Viirre et al., 1986; Chen-Huang and McCrea, 1999), (ii) cortical computations (Snyder and King, 1992), (iii) parametric changes (Green, 2000), (iv) multiplication of vestibular and eye position signals (Zhou et al., 2007) or (v) parallel linearnonlinear pathways (Lasker et al., 1999). All these models are only focused on the slow phases of VOR nystagmus.

In our recent work (Ranjbaran and Galiana, 2013a), we presented a nonlinear bilateral model for AVOR slow phases in the dark. The model is developed based on known realistic physiological mechanisms and anatomical connections including the semicircular canals, the VN and PH neural populations, motorneurons and eye plants (**Figure 2A**). Based on geometrical relations, we showed that combining monocular and vergence angle

(eye movements in opposite directions) information is sufficient to locate a target in space relative to the eyes. By assigning properly tuned nonlinear neural computations at the VN level, this slow phase model is capable of replicating target-distance dependent VOR responses that meet geometrical requirements. Nonlinear computation in neural responses, so-called *gain modulation*, exists in many cortical and subcortical areas (Salinas and Sejnowski, 2001). Different mechanisms are proposed to explain them, such as recurrent neural networks (Salinas and Abbott, 1996), changes in the synchrony of inputs to a neuron (Salinas and Sejnowski, 2000) or varying the level of background synaptic input (Chance et al., 2002). In this slow phase model (Ranjbaran and Galiana, 2013a), it is postulated that the sensitivity of the VN cells to vestibular signals modulates nonlinearly with eye position and vergence state, enabling auto-adjustment of the VOR to the set point of both eyes- a great improvement over the initially proposed model that only used ipsilateral monocular signals (Khojasteh and Galiana, 2009b). In addition to the near ideal AVOR gain modulation with target distance, the central premotor responses in that model are also consistent with experimental observations. The model also reproduces experimental observations of the VOR responses with simulated unilateral canal plugging, an emerging property. Due to nonlinearities in the sensors and premotor circuits, the model predicted a disconjugate VOR in the dark. However, prior explorations of the model behavior were examined only during high frequency head pulses or low amplitude sinusoidal rotations to remain in the range of feasible eye rotations. In order to have more relevance to the clinical VOR, we now examine the predicted responses to low frequency sinusoidal and large head rotations. This requires the implementation of a fast phase circuit to replicate realistic nystagmus patterns in the AVOR and compare simulations to experimental data.

Classically, the two phases of the VOR are believed to be generated by independent and parallel pathways as originally suggested by Chun and Robinson (1978). Based on this approach, the two phases function independently from each other and a switching strategy implements the timing of changes from one system to the other. Such a black box approach was also employed by other researchers to study VOR slow and fast phase interactions (Winters et al., 1984). However, more recent data demonstrate that slow and fast phases of the VOR share efference copies of eye position from PH and premotor cells in VN (Fukushima and Kaneko, 1995). The first model to include shared connections between the slow and fast circuit was proposed by Galiana (1991) where distinct dynamics for slow and fast phases are generated here through structural modulation. In other words some of the projections during slow phases alter their response characteristics during a fast phase: e.g., position vestibular pause (PVP) and eye head velocity (EHV) cells in the VN pause for ipsilaterally directed fast-phases (McFarland and Fuchs, 1992) and burster cells that are only active during fast phases or saccades, play an important role in facilitating response changes on premotor cells (Kitama et al., 1995). It should be noted that structural modulation does not refer to any change in anatomical connectivity, but rather to changes in the available set of active pathways (Galiana, 1991). Another comparable physiologically relevant model was also developed using realistic spiking neurons that replicated VOR nystagmus in the guinea pig with shared connections between the slow and fast circuits (Cartwright et al., 2003). However, both these models do not address VOR gain modulation with target distance nor vergence interactions. The hybrid model developed by Khojasteh and Galiana (2009b) considered VOR gain modulation; however, they considered a nonlinear block in a feedback loop in their slow phase model

**FIGURE 2 | (A)** Right and left eye positions (*ER*, *EL*) for an eccentric target at a location given by (*D*, θ) during head rotation about radius *r*. *I* is the interocular distance. **(B)** Bilateral model of slow phase horizontal AVOR in the dark. **(C)** Model structure for a rightward fast phase. Inactive projections and paused cells are indicated by dashed-gray lines. Long dashed-black lines are the centers that are only active during fast phase and solid black lines are shared

projections during both slow phase and rightward fast phase of the VOR. Silenced cross midline projections are not shown for simplicity. PVP cells are located on the opposite side from their physical location to better view the connections crossing the midline. See **Table 1** for projection weights that are shown on the connections. Projection weights that are not shown are assumed to be 1. Parts **(A)** and **(B)** partially adapted from Ranjbaran and Galiana (2013a).

to account for VOR gain modulation which resulted in variable VOR dynamics with context and head velocity profiles. Moreover, in their model only ipsilateral monocular signals are used to modulate VOR gain, thus it is not possible to test VOR gain modulation during simultaneous vergence goals and harmonic vestibular inputs.

In this paper, the fast phase circuit shares premotor centers with the bilateral nonlinear slow phase circuit previously presented (Ranjbaran and Galiana, 2013a) to form the VOR hybrid model. A novel feature is that for the first time, VOR and vergence interaction is included in a physiologically relevant hybrid model that can replicate experimental observations, i.e., modulation of the VOR gain in response to simultaneous variable vergence goals and vestibular inputs in the dark. A viable switching strategy is also implemented to trigger and stop VOR fast/slow phases, originally suggested by Galiana (1991). Simulation results are presented to evaluate the performance of this hybrid model under different rotation profiles and the results are compared

(1b)

with experimental observations. Such a model allows new interpretations of the underlying mechanisms in the VOR system and explains contradictory observations in experiments. Preliminary results are presented in Ranjbaran and Galiana (2013b).

The remainder of this paper is organized as follows. Materials and Methods in Sections 2.1 and 2.2 review briefly the reference coordinates and the previously developed slow phase model. Sections 2.3 and 2.4 describe the fast phase model and the nystagmus strategy. Simulation results in Section 3 are followed by discussion and concluding remarks in Section 4.

#### **2. MATERIALS AND METHODS**

In developing mathematical representations for the slow and fast phases of the VOR, we are modeling population responses of cells. Therefore, each element of the models represents the average behavior of a particular cell type rather than the response of any individual cell. Moreover, only firing modulation around a population resting rate is considered; biases due to resting rates are not included and a negative firing rate refers to a cell firing below its resting rate. Finally, we wish to represent the simplest model that can replicate general VOR characteristics, i.e., a minimalist approach. Adding more projections and loops between elements in a bilateral model will only affect the current assigned projection weights and not the general characteristics of the model.

#### **2.1. REFERENCE COORDINATES**

For each eye, zero position is defined as looking straight ahead at optical infinity; temporal deviations are considered positive and nasal deviations, negative. Conjugate and vergence eye positions are thus defined as *Econj* = 1 2 (*ER* − *EL*) and *Everg* = −(*ER* + *EL*), where *ER* and *EL* refer to the right and left eye angle, respectively (see **Figure 2A**).

#### **2.2. SLOW PHASE MODEL**

The original nonlinear model for slow phases of the AVOR (**Figure 2B**) is described in detail in Ranjbaran and Galiana (2013a). The input is head velocity, *sH*(*s*), sensed by semicircular canals. The canals are modeled as high-pass filters of head velocity, *V*(*s*) = *sTc sTc* + 1 , followed by a static nonlinearity on sensory modulation. The nonlinear block accounts for the mechanoneural transduction process that causes asymmetric changes in the firing rate on the primary afferents (Goldberg and Fernandez, 1971). The nonlinear block has asymmetric gains around zero (*knegative* = 0.4 and *kpositive* = 0.6) and limits the primary afferent output *VR*,*L* by saturation (+110 spikes/s) and cutoff levels (−90 spikes/s), appropriate for primary vestibular afferents with the 90 spikes/s resting rate. PVP and EHV cell populations in the VN are distinct in the model and receive sensory projections from the canals as well as efferent copies of eye position from PH. T-II in our model refers to type II neurons in the medial VN (Shimazu and Precht, 1966) that receive projections from the contralateral VN. We assume that contralateral VN projections arise from from PVP cells and form a feedback loop between the two sides of the VN (Keller and Precht, 1979). These commissural pathways play an important role in the dynamics of the VOR system (Galiana and Outerbridge, 1984). Premotor PVP and EHV cells project to motor neurons (MN) to drive the eye plants. The eye globe and muscles smooth the motoneural signals. This is represented mathematically as a low-pass transfer function. Therefore, the eye plants as well as neural filters in PH are modeled with first order low pass dynamics as *P*(*s*) = *kp sT* + 1 and *F*(*s*) = *kf sT* + 1 . Here, we assume that efference copies of the ipsilateral monocular eye position, *E*ˆ*R*,*L*, and of the vergence eye position, *E*ˆ *verg* , reach EHV cells and define their sensitivities (gain) to vestibular signals in a nonlinear fashion; i.e., *EHVR*,*L* = *gR*,*L* - *E*ˆ*R*,*L*, *E*ˆ *verg p*2*VR*,*L*, where *gR*,*L*{.} is the nonlinear sensitivity of EHV cells to vestibular afferents. These nonlinear computations account for the target distance related gain modulation of the VOR.

The equations for conjugate and vergence angles in the model are

$$E\_{conj} = \frac{k\_{\mathcal{P}}(c-1)(\mathbf{g}\_L\{\mathbf{.}\}p\_2V\_L - \mathbf{g}\mathbf{g}\{\mathbf{.}\}p\_2V\_R) - ak\_{\mathcal{P}}p\_1(V\_L - V\_R)}{2\left((c-1)(T\mathbf{s}+1) + adk\_f\right)}\tag{1a}$$

$$E\_{verg} = \frac{k\_{\mathcal{P}}(c+1)(\mathbf{g}\{\mathbf{.}\}p\_2V\_L + \mathbf{g}\mathbf{g}\{\mathbf{.}\}p\_2V\_R) - ak\_{\mathcal{P}}p\_1(V\_L + V\_R)}{(c+1)(T\mathbf{s}+1) - adk\_f}$$

Modulation of *gR*(.) and *gL*(.) changes the context gain but not the system dynamics (poles). In simpler terms, the sensory signals from the semi-circular canals are smoothed and tuned by the brainstem circuits to generate the eye movement, described by nonlinear low pass dynamics. For complete justification of the model elements and connections see (Ranjbaran and Galiana, 2013a).

The slow phase model is originally designed to replicate VOR responses in the dark with no visual cue. In order to evaluate the effect of far vs. near target flashes during sinusoidal rotation in the dark, additional inputs to trigger vergence eye movements are required. It is postulated that viewing a flashed target causes signals to be relayed to the neural filters in the PH from any cortical or brainstem center coding visuomotor error commands, such as superior colliculus (SC) (Cova and Galiana, 1996; Green, 2000). We define these visual error signals (*VeR* and *VeL*) as additional input signals to the PH (**Figure 2B**). In the absence of head velocity input, i.e., *VR*,*L* = 0, the conjugate and vergence response to *VeR*,*L* are obtained as

$$E\_{conj} = \frac{ak\_\text{p}k\_f dq (V\text{e}\_\text{R} - V\text{e}\_\text{L})}{2(\text{Ts}+1)\left((\text{c}-1)(\text{Ts}+1) + adk\_f\right)}\tag{2a}$$

$$E\_{\text{verg}} = \frac{-ak\_{\text{p}}k\_{\text{f}}dq(\text{Ve}\_{\text{R}} + \text{Ve}\_{\text{L}})}{(\text{Ts} + 1)(\text{c} + 1)(\text{Ts} + 1) - adk\_{\text{f}}} \tag{2b}$$

Assigning identical visuomotor error commands, i.e., *VeR* = *VeL*, results in pure vergence with no conjugate response due to the bilateral structure of the model. It should be noted that we are not including light conditions or continuously visible targets in the dark since the VOR dynamics will change as additional visual loops are added to the circuit (Green, 2000). This is a question for future studies.

Two sets of parameters are provided in **Table 1** that simulate two different time constants for the conjugate slow phase system obtained from Equation (1A), i.e., *Tconj* = *T*(*c*−1) *adkf* + *c* − 1 equals

**Table 1 | Numerical values of the model parameters.**

5 s or 1.2 s and *Tverg* = *T*(*c* + 1) −*adkf* + *c* + 1 from Equation (1B) equals 0.4 s with both parameter sets. By simply changing projection weights or filter gains, e.g., *c* and *d* or *kf* , new time constants for the model can be obtained. However this then requires retuning of the nonlinear surfaces at the EHVs. The canal time constant is set to *Tc* = 6 s. Nonlinear surfaces assigned to the EHV cells (Ranjbaran and Galiana, 2013a) are provided in Appendix 5.1. Here interocular distance is: *I* = 6 cm and the axis of rotation is: *r* = *rhead* = 8.8 cm.

#### **2.3. FAST PHASE MODEL**

The model structure for a rightward fast phase circuit is shown in **Figure 2C**. A leftward fast phase is generated with a mirror image of this model. Similar to the slow phase circuit, only modulations in cell populations are provided. Summing junctions are linear except for the nonlinear EHV cells (Ranjbaran and Galiana, 2013a). The bilateral structure of the slow phase system with reciprocal signals across the midline switches to a unilateral structure during fast phases. This is the result of silenced VN cells such as the ipsilateral PVPs and EHVs as well as the cross-midline projections during a fast phase (gray dashed lines and circles in **Figure 1B**, cross midline projections are not shown for simplicity). Omnipause neurons (OPN) and burster-driving neurons (BDN) as well as excitatory and inhibitory burst neurons (EBN and IBN) are included in the fast phase circuit (long dashed black lines and circles in **Figure 1**). OPNs are located near the midline of the pons and act as triggers for the initiation of fast eye movements in all directions (Scudder et al., 2002). OPNs discharge at high firing rate and exert a tonic inhibition on premotor BNs during slow eye movements and fixation. Prior to a fast eye movement or saccade, OPNs cease firing, remain silent during the saccade, and resume firing as the saccade ends (Yoshida et al., 1999). OPNs receive projections from the SC as well as projections from cells in the medial VN (Ito et al., 1986). In our model, it is assumed that the projections from the medial VN, specifically PVPs, to the OPNs play a role in triggering and ending the fast phases (see Section 2.4). BDNs are located below the PH and are found to be excited by contralateral horizontal head rotation and they send projections to contralateral BNs (Kitama et al., 1995). We assume BDN excitation is a result of an excitatory vestibular drive that comes from contralateral vestibular-only (VO) cells. BDNs also modulate with a PH response (eye position efference copies) to close the loop and shape bursts during fast phases (Kitama et al., 1995). Therefore, the output signal from left BDNs during rightward fast phase is: *BDNL* = *m* × *VR* − α × *E*ˆ*R* where α is the projection weight from the PH to contralateral BDN (connections are simplified in **Figure 1**). BDNs project to contralateral EBNs and IBNs that are located in the reticular formation and have monosynaptic connections to abducens motoneurons. EBNs send excitatory projections to ipsilateral MNs while IBNs with similar firing patterns send inhibitory projections to contralateral MNs. In our model, the same projections to MNs are sent to PH neurons that produce efferent copies of eye position for VN cells (Fukushima and Kaneko, 1995).

In order to achieve faster dynamics in the fast phase circuit, it is assumed that the feedback loops including PVP and EHV cells between the VN and PH nuclei change their net sensitivity direction as originally suggested in Galiana (1991). This can be a result of competition between parallel inhibitory and excitatory projections from PVP and EHV cells whose balance is modified by EBN/IBN effects. Projections from burst neurons to VN are studied in Igusa et al. (1980) and have been used in former models and studies (Curthoys, 2002; Cartwright et al., 2003). We assume that the population response of the VN cells is a mixed combination of individual inhibitory and excitatory projections. During slow phases, the population response is dominated by the excitatory PVPs and inhibitory EHVs as the burst neurons are silenced. During a fast phase, however, as the burst neurons become dis-inhibited, along with silencing of the ipsilateral PVPs and EHVs, the inhibitory projections dominate the response of PVPs and excitatory projections dominate the response of EHVs; thus, the net projections from contralateral PVPs and EHV populations appear to change their direction of sensitivity. Similar incrementingdecrementing behavior is also observed on the BDN activity during slow-fast intervals (Ohki et al., 1988). The change in effective connectivity assumed here does not conflict with basic knowledge about the neural firing patterns. Thus, in our model during a fast phase, PVPs and EHVs contralateral to the fast phase direction change their sensitivity and their activity profiles decay.

As in the slow phase model, the eye plants and neural filters in the PH remain as first-order low-pass dynamics as *Pf*(*s*) = *kpf Ts* + 1 and *Ff*(*s*) = *kff Ts* + 1 . EHV cells are silenced during ipsilateral fast phases and are active during contralateral fast phases with the same nonlinear sensitivity to the ipsilateral canal signal, i.e., *g*(.) (Ranjbaran and Galiana, 2013a). Due to the unilateral structure of the fast phase of the VOR, distinct monocular dynamics for eye responses are obtained during fast phases, i.e., during a rightward fast phase (see Appendix 5.2)

$$E\_R = \frac{k\_{\mathcal{pf}}(-ap\_1V\_L + mb\_EV\_R)}{T\text{s} + 1 + k\_{\mathcal{ff}}(ad + b\_E\alpha)}\tag{3a}$$

$$E\_L = \frac{k\_{\rm{pf}} \left(\beta\_1 V\_L + \beta\_2 V\_R\right)}{(1 + T\text{s})(T\text{s} + 1 + k\_{\rm{fl}}(ad + b\_E\alpha))}\tag{3b}$$

where

$$\begin{cases} \beta\_1 = (1 + \text{Ts} + k\_{\overline{\mathcal{H}}}(ad + b\_E \alpha)) \mathbf{g} \{ .\} p\_2 - b\_I \alpha k\_{\overline{\mathcal{H}}} a p\_1 \\ \beta\_2 = -(1 + \text{Ts} + k\_{\overline{\mathcal{H}}}(ad + b\_E \alpha)) m b\_I + b\_I \alpha k\_{\overline{\mathcal{H}}} m b\_E \end{cases} (4)$$

During a leftward fast phase, the dynamic equations for *ER* and *EL* are obtained by switching the *R* and *L* subscripts in Equations (3A,B). The above equations imply that monocular eye trajectories have different dynamics during rightward and leftward fast eye movements.

The model parameters (**Table 1**) are selected to preserve the stability of the fast phase system with a small time constant.

#### **2.4. STRATEGY FOR NYSTAGMUS**

So far, the slow and fast phase models with shared connections are described. However, an important feature of the VOR is the switching mechanism between these two phases. The linear range of the VOR is improved by nystagmus as eye excursions are kept inside a reasonable limit. Therefore, a proposed switching strategy is based on limiting eye deviations by avoiding cut-off and saturation limits in the responses of premotor neurons (Galiana, 1991). It is known that the activity of OPNs acts as a logical circuit to trigger and end a fast phase. In this model, the OPN circuit constantly monitors the output of PVPs. If the firing rate of PVPs on one side reaches a threshold (ON-Th spikes/s), a fast phase is triggered ipsilateral to the PVPs' side. During the fast phase, the ipsilateral PVPs and EHVs are silenced and the contralateral PVPs and EHVs change direction and decay as explained in the fast phase circuit structural modulation. The fast phase ends as the firing rate of the contralateral PVPs decays to a second threshold (OFF-Th spikes/s). The fast phase intervals are therefore generated in the same direction as head movement and eye position signals are kept below their physical limits. ON-Th and OFF-Th control the frequency of fast phases and their duration. For instance, with fixed model parameters and dynamics, increasing ON-Th results in later triggering of fast phases and lowering the OFF-Th leads to longer fast intervals. We have imposed a refractory period of 20 ms in the model after switching back to a slow phase to enforce a minimum time interval before triggering a new fast phase.

The performance of the model under different conditions are provided next. All simulations were performed using MATLAB Simulink (The MathWorks Inc., USA), with a first order Euler approximation and a step size of 1 ms.

### **3. RESULTS**

The model is designed to simulate the human AVOR responses during yaw rotations around a vertical axis centered on the head. We focus here on the global, behavioral aspects of the AVOR model rather than on individual components. PVP and EHV firing behavior is previously addressed in Ranjbaran and Galiana (2013a).

#### **3.1. RESPONSE TO SINUSOIDAL ROTATION IN DARKNESS**

**Figure 3** depicts the response of the hybrid model with *Tconj* = 5 s at two different rotation frequencies: 1/6 Hz (A,B) and 1/2 Hz (C,D) with velocity peaks of 180 degree/s. As in experimental observations, the number of fast phases *per cycle* decreases for higher frequency sinusoidal head rotations. In other words, fast phases are triggered more often during low frequency head rotations. This is due to the band pass characteristics of the central

neurons in the VOR pathway. At lower frequencies the gain of central neurons is higher which increases the possibility of exceeding their firing thresholds and triggering a fast phase. Furthermore, at a given rotation frequency, fast phases are more frequent at the higher head velocity levels; consistent with experimental observations (Buettner et al., 1978).

We have also compared our model performance with *Tconj* = 1.2 s in response to a specific rotation profile (180 degree/s at 1/6 Hz) where binocular records are available in our archive; for details of the experiment see (Khojasteh and Galiana, 2009a). **Figure 4** depicts the recorded conjugate eye position (A) and eye velocity (B) (gray), as well as our model responses to this rotation stimulus (black). Clearly, general nystagmus characteristics in the simulation and data are similar: the amplitude of conjugate eye position, the number of fast phases and their timing, suggest that the switching mechanism is plausible.

Given nonlinear canals and nonlinear premotor (EHV) computations, eye movements are disconjugate and a vergence component is now present in the response of our hybrid model to head perturbations. This vergence response shows a carrier frequency that is twice that of the stimulus, as also seen in the experimental data (**Figure 4C**). The peak-to-peak amplitude of the vergence component is greater than EOG resolution limits, suggesting that it cannot be a result of inappropriate calibration in the binocular recording. Instead, as predicted from the nonlinear model, this vergence component can be a direct result of nonlinearities at the premotor and sensory levels. It should be noted that here we did not attempt to identify a system directly from eye recordings but rather compare the general characteristics of recorded AVOR and our model response.

#### **3.2. CONTEXT DEPENDENT RESPONSE TO SINUSOIDAL ROTATION**

The experimental work of Paige et al. (1998) studied the role of fixation distance in adjusting the gain of the VOR during *sinusoidal* angular head rotation. Here, we will test the response of our model during sinusoidal rotations while fixating on a near or far flashed target. Vergence movement as a result of fixating a target in the model is obtained by assigning proper visuomotor error commands, i.e., *VeR*,*L*.

Starting with zero initial conditions (i.e., looking straight ahead at optical infinity), *ER*(0) = *EL*(0) = 0, *VeR* = *VeL* are set to replicate a flashed target in the dark appearing in the sagittal plane between the eyes. This flashed target appears 5 s after the start of head rotation, at *D* = 43 cm from the eyes requiring 8 degree of vergence (given the interocular distance of *I* = 6 cm; see Ranjbaran and Galiana, 2013a for geometrical relations). At *t* = 10 s, a new flashed target appears at *D* = 21 cm that requires 16 degree vergence. At *t* = 15 s a far flashed target appears to reset vergence to zero.

The simulation results are obtained with model parameters where *Tconj* = 1.2 s and *Tverg* = 0.4 s (see **Table 1**). In the absence of head rotation, there is no conjugate eye movement, only a vergence response (**Figure 5**). During sinusoidal head rotation at 0.5 Hz with 120 degree/s peak velocity and the same visuomotor inputs, the conjugate and vergence AVOR responses are appropriate (**Figure 6**). The gain of the VOR (peak of envelope of eye velocity/peak head velocity) is near unity during the first 5 s with no visuomotor response as expected. During *t* = 5 → 10 s interval, the first flashed target, causes a VOR gain increase from unity to compensate for the visuomotor command and the vergence movement. As the second target appears at *t* = 10 s, a still larger vergence is required and therefore the VOR gain also increases further. As the final far target appears at *t* = 15 s, the vergence position decays to zero (looking far ahead) and the AVOR gain decreases smoothly to default unity. This is in agreement with the experimental observations (Paige et al., 1998; Viirre et al., 1986). Note that the pure vergence response seen in **Figure 5** is combined with vergence modulations during each fast phase of nystagmus in **Figure 6**: this has been seen by Sylvestre et al. (2002) during visual disjunctive saccades, an emerging property.

We also tested the effect of head rotation frequency on the gain of the VOR while fixating central flashed targets at different depths. According to observations by Paige et al. (1998), the AVOR gain increases with rotation frequencies while fixating an *imaginary* earth fixed target in darkness. Moreover, in plots of the resulting VOR against concurrent vergence, they report that both the slope, and the intercept of this graph at 0 vergence, increase with rotation frequency (see Figure 8 in Paige et al., 1998). We emphasize the context of imaginary target, since our model does

eye position (degree) and visuomotor command (*VeR* + *VeL*, gray). Vertical dashed lines mark the time of change in the visuomotor command.

not include vision related loops in the light or constant visual input, but rather a vergence cue given by a flashed target.

To replicate this experiment, the hybrid model with *Tconj* = 1.2 s is simulated with sinusoidal head rotation in the range of 1/6–4 Hz, with 30 degree/s peak head velocity while fixating central flashed targets at one of these distances: *D* = [2000 85.9 42.9 28.54 21.34 17.01] cm. Given the interocular distance and the rotation radius, the vergence angles for these target distances are: [0 4 8 12 16 20] degree, respectively. **Figure 7** describes the effects of rotation frequency. The AVOR gain during low frequency rotation, 0.5 Hz, is closer to the ideal gains obtained from geometrical relations (for detail see Ranjbaran and Galiana, 2013a) compared to higher frequency rotation at 4 Hz (**Figure 7A**). Moreover, the slope and intercept plots (**Figures 7B,C**) show an increase with increasing frequency. These emerging results are in agreement with experimental observations of Paige et al. (1998).

#### **3.3. DO CHANGES IN VOR ANTICIPATE CHANGES IN VERGENCE ANGLE?**

The work of Snyder and King (1992)investigated the contribution of the vergence angle to VOR performance. In their experiments, they measured eye velocity in rotating monkeys while the vergence angle was required to change by flashing targets at different distances. Their results demonstrated that the VOR gain changed toward its correct value for the new target distance, before the correct vergence was acquired. They concluded that VOR gain modulation by target distance anticipates changes in vergence angle; thus, the binocular vergence angle alone, derived either from proprioception or from efference copy of motor command, is not sufficient to drive VOR modulation. They suggested that the transient discharge of gaze velocity Purkinje cells in the flocculus, associated with changes in vergence angle early enough, could drive VOR gain modulation with target distance (Snyder and King, 1992).

Here, we replicate their experiment with our nonlinear model. The visuomotor inputs, *VeR*,*L* are applied after *t* = 100 ms such that a vergence movement is generated from 0 degree to 8 degree. This emulates a subject initially fixating on a far target straight ahead, and then fixating on a near central target at *D* = 43 cm. During this vergence movement lasting ≈ 2.5 s in our model, a vestibular input, i.e., a pulse of head velocity (accelerating with 500 degree/s2 to 30 degree/s, maintained for 40 ms and then decelerated) is added to the model and the peak eye velocity response is measured. As done by Snyder and King (1992), this vestibular input is applied at different times from 0 to 2.5 s in steps of 50 ms or 100 ms during the time course of the vergence movement (one pulse in each trace). Given the brief and small head perturbation, no fast phase is triggered. **Figure 8** shows that the velocity of the eye movement evoked by VOR changed smoothly over the course of the 8 degree convergence. Both VOR peak eye velocity and vergence angle were normalized and replotted as functions of time to resemble the experimental results obtained by Snyder and King (1992) (their **Figure 3**). Our model simulations replicate their observations: VOR gain changes lead the vergence angle changes, and the change in VOR was completed before the vergence angle reached its goal.

Contrary to the suggestion by Snyder and King (1992), here, this observation is a result of VOR gain modulation by the *efference copies* of the vergence and monocular angles at the premotor level. In our nonlinear model, the effect of the required change in vergence appears in the efference copies *E*ˆ*R*,*L* before *ER*,*L* via visuomotor projections to the PH. Consequently, the projections to the EHVs from the PH carry early vergence information to modulate the VOR gain, before the vergence movement actually begins or completes. It may appear that VOR changes anticipate vergence angle changes, but instead we suggest that the efference copies of the vergence angle modulating the VOR gain carry the vergence command before vergence responses appear at the behavioral level.

#### **3.4. RESPONSE TO STEPS IN HEAD VELOCITY**

Raphan et al. (1979) explored extensively the characteristics of VOR nystagmus during and after steps of passive head velocity,

both in the dark and in the light. Eye velocity profiles in man and monkey decay to zero in the dark if the rotation interval exceeds 2–3 times the time constant of the canals. This is expected from a high-pass sensor. The simulations below focus on nystagmus in the dark, the context of our current model, using steps of head velocity of variable amplitude and duration.

#### *3.4.1. Per and post rotatory nystagmus*

The simulations in **Figure 9** replicate the main characteristics of VOR nystagmus during and after steps of head velocity. First, for long 45 s steps of rotation, it is clear that the post-rotatory nystagmus velocity appears equal in magnitude but opposite in direction to that during the rotation. In addition, the nystagmus velocity peak scales in both per- and post-rotation with the amplitude of the head velocity (**Figures 9B,C**). Second, for short duration rotations (**Figures 9D,E**), the initial post-rotatory nystagmus in the opposite direction is reduced in magnitude as the interval of rotation shortens. This is expected since the high-pass dynamics of the sensor will cause the *change* of nystagmus velocity to remain constant, if measured from the current eye velocity at the moment rotation ceases (see arrow in **Figure 9D,E**). As found in experimental data (Raphan et al., 1979), the frequency of fast phases in all cases is also modulated by the concurrent level of slow phase eye velocity.

#### *3.4.2. Dynamics of nystagmus decay*

Raphan et al. (1979) also studied the decay rate of nystagmus velocity. As commonly done in the literature, they evaluated the dynamics of the VOR slow phase velocity by removing fast phases and replacing the gaps by interpolation: the reconstructed *envelope* was deemed to represent VOR dynamics, fitted with exponentials. The main result is that nystagmus decay appears much slower than the underlying slow-phase system (≈ 15 s vs. 4–6 s canal), hence the term *velocity storage* in the VOR coined by Raphan et al. However, an envelope fit ignores the contribution of initial conditions introduced at the start of each slow phase segment, biasing estimates of the slow phase time constant. To illustrate, **Figure 10** provides the hybrid model response to a step of −250 degree/s in head velocity, with a canal time constant of 6 s and a conjugate slow-phase time constant of 1.2 s (see **Table 1**). The slow phase central time constant is intentionally low to highlight the effects, but these hold whenever there is nystagmus, especially at very low frequencies like steps. In **Figure 10A**, the *envelope* of slow phase velocities decays with a time constant of 5.55 s, despite slow phase central dynamics of 1.2 s. Such a response is often seen in unilateral vestibular patients. As discussed for sinusoidal rotations in Galiana (1991), ignoring the effect of nystagmus results in biasing the estimated conjugate VOR dynamics. There is a plateau-like response in the initial nystagmus velocity also seen by Raphan et al. (1979) at higher head speeds. Here it is caused by nonlinearities in the canal sensitivity, now exceeded by the input range. In addition, with the hybrid model, we predict the appearance of vergence nystagmus (**Figure 10B**) during step rotations in the dark. In order to extend *velocity storage* beyond both canal and central time constants, it is sufficient to incorporate the resting rates of sensors and central components (Galiana, 1991); at this time we only include modulations at all sites about resting rates, so only the central time constants can be *masked* during nystagmus.

## **4. DISCUSSION**

This paper introduces a hybrid nonlinear model to replicate human AVOR nystagmus in the dark. This bilateral model includes nonlinear sensors as well as nonlinear surfaces assigned to EHV cells to account for the target distance dependence of the VOR. It is shown that vergence can appear with both vestibular and visual depth stimuli. A physiologically relevant fast phase circuit and a nystagmus strategy are imbedded to generate nystagmus automatically and extend the functional range of the AVOR. In our former work (Ranjbaran and Galiana, 2013a), a comprehensive study was performed on the slow phase aspects of the VOR and their characteristics under different conditions. Here, we evaluated the performance of the hybrid model through simulations in response to passive head rotations with or without flashed visual goals. For the first time, a hybrid model based on known brainstem connections replicates different VOR characteristics, consistent with experimental observations. The comparison between the hybrid model simulation responses and experimental observations are qualitative at this stage, showing the capacity to explain reported behavior. Clearly experiments on individual subjects results in distinct numerical responses which require retuning of the parameters in the hybrid model.

**to constant head velocity rotation. (A)** Conjugate position (degree) and **(B)** Conjugate vel. (degree/s) for for 45 s head rotation at −60 degree/s. **(C)** same as **B**, but doubling the speed of head rotation to −120

now with short intervals of 10 and 3 s respectively. **(A)** black: conjugate position and gray: negative head vel./3 . **(B–E)** black: conjugate vel., gray: negative head vel.

Simulated nystagmus patterns replicate reported experimental observations (**Figure 3**) and trajectories that resemble human data (**Figure 4**). This suggests that the switching mechanism in our model is both plausible and testable with lesions and new inputs.

### **4.1. DISCONJUGACY OF THE AVOR:**

Contrary to common belief, the AVOR is not purely conjugate in the dark; binocular recordings during sinusoidal rotations in darkness confirmed a vergence component in the AVOR (Khojasteh and Galiana, 2009a). In our model, this vergence component is a result of nonlinear sensors as well as nonlinear premotor cell responses that account for context dependent VOR responses. This suggests that local nonlinearities in the VOR circuit are the underlying mechanism for the disconjugate VOR in the dark.

## **4.2. VESTIBULAR VERGENCE INTERACTIONS:**

In addition to the vestibular input, i.e., head movement, visuomotor commands are included to enable vergence movements in response to *flashed* targets in the dark. Simulations show the effect of vergence goals during sinusoidal rotations: they confirm that the context dependency of the AVOR gain in the model is preserved with nystagmus and variable vergence goals. AVOR gain dependency on rotation frequency is also in agreement with experimental observations (Paige et al., 1998), an emerging property.

It appears experimentally that AVOR gain modulation with target distance precedes changes in vergence (Snyder and King, 1992), which questions vergence itself as the drive for AVOR gain modulation. We replicated this experiment (Snyder and King, 1992) with our model and found the same result: AVOR gain changes anticipate or precede the vergence profile

(**Figure 8**). We conclude that AVOR gain modulation using efference copies of the vergence angle can support this anticipatory VOR modulation, even in the dark: we postulate that efference copies with visuomotor inputs affect EHV cells immediately, and so modulate the AVOR gain before the behavioral vergence is fully executed.

### **4.3. AVOR DYNAMICS DURING STEPS OF HEAD VELOCITY:**

Per and post rotatory nystagmus characteristics are influenced by the peak velocity and duration of the stimuli. The proposed model replicates these data patterns (**Figure 9**). Given the switching aspect of nystagmus, we demonstrate that *envelope* measures provide biased estimates of slow-phase dynamics (**Figure 10**). So an important goal is to develop algorithms that provide unbiased estimates of nystagmus dynamics.

### **4.4. TESTABLE PREDICTIONS:**

The goal of modeling a sensory-motor system is revealing potential strategies in the brain to control motion and to gain insight for clinical applications. Since the modeling results are fully consistent with available experimental data, the model structure warrants further study. Several assumptions or predictions remain to be verified: (i) Assumptions regarding anatomy of the VOR:

• The presence of premotor (e.g., PVP) projections to OPN cells, to enable the proposed switching strategy.

(ii) Predictions regarding VOR dynamics and behavior:

In summary, we explored AVOR online gain modulation with target distance by introducing a physiologically relevant hybrid nonlinear model. We proposed local nonlinear computations at VN levels to account for the gain modulation of the VOR with context. It is likely that this hypothesis could also support long term adaptation or lesion compensation in the VOR. Furthermore, this hybrid model, given the realistic aspect of its simulated data, can also be used to generate virtual data for validation of algorithms that classify nystagmus segments and identify reflex dynamics.

# **FUNDING**

This work has been supported by Canadian Institutes of Health Research (CIHR), Natural Sciences and Engineering Research Council of Canada (NSERC) and Fonds de recherche du QuÉbec (FQRNT).

# **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 14 July 2014; accepted: 14 January 2015; published online: 09 February 2015. Citation: Ranjbaran M and Galiana HL (2015) Hybrid model of the context dependent vestibulo-ocular reflex: implications for vergence-version interactions. Front. Comput. Neurosci. 9:6. doi: 10.3389/fncom.2015.00006*

*This article was submitted to the journal Frontiers in Computational Neuroscience. Copyright © 2015 Ranjbaran and Galiana. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **APPENDIX**

#### **5.1. NONLINEAR SURFACE**

The nonlinear surface assigned to the EHV cells, *gR*,*L*, is obtained as a 3th order polynomial of *x* = *E*ˆ*R* for the right EHV and *x* = *E*ˆ*L* for the left EHV and a 1st order polynomial of *y* = *E*ˆ *verg* according to the following equation:

$$\mathbf{g}\_{R,L} = m\_0 + m\_1 \mathbf{x} + m\_2 \mathbf{y} + m\_3 \mathbf{x}^2 + m\_4 \mathbf{x} \mathbf{y} + m\_5 \mathbf{x}^3 + m\_6 \mathbf{x}^2 \mathbf{y} \tag{A1}$$

For the model with *Tconj* = 1.2 s, the coefficients in Equation (A1) are: *m*0 = 2.59, *m*1 = −8.051*e* − 5, *m*2 = 0.12, *m*3 = −4.8*e* − 6, *m*4 = 1.52*e* − 5, *m*5 = −4.12*e* − 6, *m*6 = −1.19*e* − 7.

For the model with *Tconj* = 5 s, these coefficients are re-tuned to: *m*0 = 1.68, *m*1 = −6.43*e* − 5, *m*2 = 0.09, *m*3 = −3.84*e* − 6, *m*4 = −1.21*e* − 5, *m*5 = −3.29*e* − 6, *m*6 = 9.54*e* − 8.

#### **5.2. FAST PHASE MODEL EQUATIONS**

In this section, the equations to obtain the fast phase dynamic relations in Equations (3A,B) are described. Subscripts *R* and *L* refer to the right and left side of the brainstem, respectively. PVP, EHV ,BDN, EBN, and IBN here refer to the net output from these cell populations. The parameters *a*, *p*1, *p*2, *d*, *m*, *bI*, *bE*, and α define weight of the projections between cell types or brainstem centers according to **Figure 2C**. Lower case letter *s* is the complex Laplace variable.

During a rightward fast phase, ipsilateral *PVP*s and *EHV*s and contralateral *IBN*s and *EBN*s are silenced. Signals from *BDNL* = *m* × *VR* − α × *E*ˆ*R* are projected to *EBNR* and *IBNR*; therefore

$$\begin{cases} EBN\_R = mV\_R - \alpha \triangle\_R & EBN\_L = 0\\ IBN\_R = mV\_R - \alpha \triangle\_R & IBN\_L = 0 \end{cases} \tag{A2}$$

PVPs receive projections through linear summation of the ipsilateral canal *V*, contralateral PVPs and the PH. EHVs receive ipsilateral canal projections *V* as well as efference copies of ipsilateral eye position and vergence angle. The net output of EHVs ,however, is scaled by a nonlinear gain that modulates the weight of canal projections according to the concurrent ocular angles (Ranjbaran and Galiana, 2013a); therefore,

$$\begin{cases} \text{PVP}\_{L} = d\triangle\_{R} + p\_{1}V\_{L} & \text{PVP}\_{R} = 0\\ \text{EHV}\_{L} = \text{g}\_{L}p\_{2}V\_{L} & \text{EHV}\_{R} = 0 \end{cases} \tag{A3}$$

Signals from VN and BN are added at MNs and then relayed to the eye plant *Pf*(*s*) = *kpf Ts* + 1 to generate ocular movements:

$$\begin{cases} E\_R = (-aPVP\_L + b\_E EBN\_R) \left(\frac{k\_{pf}}{Ts + 1}\right) \\ E\_L = (EHV\_L - b\_IIBN\_R) \left(\frac{k\_{pf}}{Ts + 1}\right) \end{cases} \tag{A4}$$

Eye position efference copies *E*ˆ*R*,*L* are available through projections from PH with similar dynamics as the eye plant, *Ff*(*s*) = *kff Ts* + 1 ; therefore,

$$
\triangle\_{R,L} = \frac{k\_{\mathcal{ff}}}{k\_{\mathcal{pf}}} E\_{R,L} \tag{A5}
$$

By substituting Equation (A5) into Equations (A2) and (A3) and combining with Equation (A4), one can obtain the dynamic equations provided in Equations (3A,B), to describe *ER*,*L*.

# Major remaining gaps in models of sensorimotor systems

#### Gerald E. Loeb1 \* and George A. Tsianos 2

1 Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, USA, 2 L-3 Applied Technologies Inc., HEM Division, San Diego, CA, USA

Experimental descriptions of the anatomy and physiology of individual components of sensorimotor systems have revealed substantial complexity, making it difficult to intuit how complete systems might work. This has led to increasing efforts to develop and employ mathematical models to study the emergent properties of such systems. Conversely, the development of such models tends to reveal shortcomings in the experimental database upon which models must be constructed and validated. In both cases models are most useful when they point up discrepancies between what we think we know and possibilities that we may have overlooked. This overview considers those components of complete sensorimotor systems that currently appear to be potentially important but poorly understood. These are generally omitted completely from modeled systems or buried in implicit assumptions that underlie the design of the model.

#### Edited by: