Vincent C. K. Cheung, The Chinese University of Hong Kong, China

#### Reviewed by:

William Zev Rymer, Rehabilitation Institute of Chicago, USA John Kalaska, Université de Montréal, Canada

#### \*Correspondence:

Gerald E. Loeb, Department of Biomedical Engineering, University of Southern California, Denney Research Center (DRB)/Downey Way, Los Angeles, CA 90089, USA<EMAIL_ADDRESS>
> Received: 16 March 2015 Accepted: 21 May 2015 Published: 04 June 2015

#### Citation:

Loeb GE and Tsianos GA (2015) Major remaining gaps in models of sensorimotor systems. Front. Comput. Neurosci. 9:70. doi: 10.3389/fncom.2015.00070 Keywords: sensorimotor control, sensorimotor learning, sensorimotor integration, sensorimotor systems modeling, biological neural networks

# Introduction

"When you can measure what you are speaking about, and express it in numbers, you know something about it, when you cannot express it in numbers, your knowledge is of a meager and unsatisfactory kind; it may be the beginning of knowledge, but you have scarcely in your thoughts advanced to the stage of science."

—William Thomson, 1st Baron Kelvin, 1883

Quantitative methodology has gradually replaced qualitative "butterfly collecting" in biology. Lord Kelvin's comment was motivated primarily by its importance for reductionistic science, whereby hypotheses about deep and unobservable structure and function can be tested according to their often subtle effects on measureable phenomena. More recently, the challenge of science has more often been too much rather than too little quantitative data, and too many well-understood but complex mechanisms whose interactions defy intuitive understanding of how complete systems actually function. This problem combined with the rapid advance of computing power has led to rapidly increasing interest in systems modeling. For the 21st century, one might replace "measure" with "model" to update Lord Kelvin's 19th century exhortation.

Modeling is a never-ending task. It is most useful when it reveals discrepancies between what we think we know about a component or a complete system and how it actually behaves experimentally. The existence of a model (whether explicit or implicit) is an inspiration to experimentalists to identify phenomena or conditions for which the model predictions are in error—the basic scientific method of hypothesis formation and falsification. Assuming the experimental results are valid, the discrepancy can only be resolved by correcting and usually complexifying the model to reflect better the properties of the physical system being modeled. As more experimental data become available, it becomes possible, indeed necessary, to break complex systems into subsystems that can be studied and modeled in relative isolation. This leads to a proliferation of models of subsystems that must eventually be combined but that may be at very different stages of development and accuracy. This problem is amplified by the natural tendency of scientists to focus on and further refine those subsystems that have already yielded to their efforts. Put simply, we have a lot of quantitative information incorporated into accurate models of some subsystems but little or none about other subsystems that are likely to be just as important to overall function or about most of the connections between subsystems.

The nature of the modeling challenge depends on where a subsystem is located in a system that is inherently hierarchical (**Figure 1**). The subsystems in the spinal cord and musculoskeletal system have mostly been studied in isolated or highly reduced preparations. Such subsystems are amenable to "bottom up" modeling strategies in which individual elements are characterized and then combined into larger models. Shortcomings in the constituent models tend to arise because it is difficult for the experimenter to observe or create the full range of natural conditions of use. The subsystems in the brain have been studied mostly in intact, naturally behaving animals. Those subsystems are generally modeled using "top-down" strategies in which the form of the model is intuited from observable phenomena of the whole system. Shortcomings in those models tend to arise because the experimenter must make assumptions about what is or is not happening in the other subsystems that are present. This overview is not an encyclopedic review of what models already exist. Rather it attempts to prioritize those subsystems that currently appear to be both important and relatively poorly modeled and to identify opportunities to rectify this.

# Systems Model Architecture

"It can scarcely be denied that the supreme goal of all theory is to make the irreducible basic elements as simple and as few as possible without having to surrender the adequate representation of a single datum of experience." (Frequently paraphrased as "A scientific theory should be as simple as possible, but no simpler.") —Albert Einstein, 1933

Biological systems (and their models) tend to be much more inherently complex than the physical systems that Einstein had in mind. This is an inevitable consequence of their gradual evolution over hundreds of millions of years of intense competition with other species. Relatively simple mechanisms such as actin-myosin binding for force generation and a stretch reflex to stabilize posture can have their essence captured by mathematical curve-fitting (Hill's equation) and engineering metaphor (servo-control), respectively. Real organisms, however, achieve their competitive performance by huge and disorderly elaborations of those underlying mechanisms. This poses a potential conflict of interest between the modeler, who often aspires to simple and elegant models, and the experimentalist, who needs to capture realistic performance.

The modeler must decide what performance is to be modeled, thereby defining the functional elements that must be included in the model system. This is itself a form of modeling that is fraught with opportunities for errors of omission. **Figure 1** provides one framework for the restricted set of learned, voluntary sensorimotor behaviors such as manipulating objects. Smaller and simpler models have been used with substantial success to account for preprogrammed behaviors such as locomotion, breathing, and mastication. Larger and more complex models will be required to account for multimodal behaviors such as eye-hand coordination.

The framework in **Figure 1** can be divided hierarchically and phylogenetically into brain, spinal cord and peripheral/mechanical subsystems. The first organisms to achieve motility evolved by steadily enhancing their abilities to make movements that were environmentally responsive, mechanically stable, and energetically efficient. Sensory transduction and electromechanical activation were already well on their evolutionary paths before there were recognizable nervous systems at all. As organisms became larger and more complex mechanically, coordination required centralization of sensorimotor connections in invertebrate ganglia that eventually coalesced into the vertebrate spinal cord. Systemslevel modeling is often applied to behaviors that are learned by the brain but implemented by the "lower" but highly evolved subsystems, whose intrinsic properties thereby define the control problems that the brain must solve. Model systems that omit or substantially simplify components are implicitly hypothesizing that those components do not make a significant contribution to the observed behavior, even though they are known to be necessary and perhaps even sufficient for many other behaviors of the same organisms. Unfortunately, such omissions more often reflect the unavailability of computational models rather than such a plausible hypothesis.

Explicit or implicit models that have been incorrectly simplified lead to progressively more complex and implausible models as they try to account for new data. This is analogous to the problem of "the music of the spheres" in which an intuitive and simple earth-centered universe requires ever more artificial structure to account for the observed motion of the planets. Such a problem may be starting to emerge in motor learning models that assume that the motor cortex is directly responsible for converting visual targets in extrapersonal space into sequences of muscle activation that cause limbs to reach to those targets. Such models often assume that the cortex learns an "internal model" of the musculoskeletal plant and inverts that model to compute the commands required to perform a given task. As Nikolai Bernstein pointed out (Bernstein, 1967) (English translation of 1934 publication in Russian), such a computational problem is ill-posed because of redundancy. There are usually more muscles and degrees of freedom than required to perform the task, so the computation requires either an arbitrary constraint on allowable strategies (d'Avella et al., 2006) or an optimization criterion such as minimizing effort that will result in a singular solution (reviewed in Loeb, 2012). If there were a single internal model in a single place subject to a single computational strategy, then one would expect this to give rise to rather consistent patterns of

sensorimotor adaptation and learning. Experiments on learning and adaptation instead reveal many different "rules" whereby subjects learn and forget how to deal with distortions in the visual space, loads on the limb, and changes in posture and muscle function and performance criteria, as well as interactions among those variables (see Shadmehr and Mussa-Ivaldi, 1994; Gandolfo et al., 1996; Krakauer et al., 2000, 2006; Baraduc and Wolpert, 2002; Mattar and Ostry, 2007, 2010; Pearson et al., 2010; Brayanov et al., 2012; Berniker et al., 2014). Substantial evidence has shown, however, that internal representations of behavior are neither intuitive nor simple (Brayanov et al., 2012) and that sensorimotor learning generalizes poorly in many situations (Gandolfo et al., 1996; Mattar and Ostry, 2007, 2010; de Rugy et al., 2012b; Coelho et al., 2013; Berniker et al., 2014), which is inconsistent with the predictions of a singular internal model. If these behaviors were actually the result of several different sensorimotor subsystems, each with their own relatively simple rules, the inferred models might actually be simpler as well as more realistic. Such a collection of interacting subsystems is known to subserve control of gaze, which consists of anatomically distinct subsystems for classes of behavior such as voluntary and involuntary saccades, smooth pursuit, and reflexive stabilization, and the coordination of eye and head movement to achieve them. The anatomical structures responsible for gaze control include the mesencephalic tectum, pontine, and other brainstem nuclei and cerebellum as well as sensory and motor cortical areas, all of which also have strong sensory and motor connections with the limbs.

# Gaps in Peripheral and Mechanical Models

Models of the musculoskeletal system include several different types of models that must be combined to generate their complete input/output properties.

## Musculoskeletal Dynamics

Control engineers divide any system into a plant and a controller and describe the tasks to be performed using cost functions that weight the relative importance of quantifiable state variables such as accuracy, time to completion and energy consumption. The job of the sensorimotor nervous system (controller) is to compute and implement control signals that cause the musculoskeletal system (plant) to achieve desired performance. The mechanical dynamics of the plant obviously constrain the set of useful solutions for the controller. More complex plants make it more difficult for engineers to develop mathematical models that can be used to compute solutions, but they do not necessarily make the plant inherently more difficult to control. For example, the complex intrinsic properties of muscles (e.g., dependency of force output on instantaneous muscle length and velocity) may provide rapid stabilizing effects that compensate in part for the relatively long delays inherent in signal transmission in a neural controller (Hogan, 1984; Loeb et al., 2002) but see (Crevecoeur and Scott, 2014). Nevertheless, the mechanical dynamics of multilinked skeletal segments such as a limb or vertebral column are inherently complex and can result in highly unstable conditions that must be prevented by a controller with widely distributed inputs and outputs (Lackner and DiZio, 2000).

Sophisticated algorithms for solving the mechanical dynamics of free-body systems with mechanical constraints have been widely applied to the biomechanics of musculoskeletal systems. These work reliably when the biological architecture lends itself to decomposition into discrete and independent entities—the inertial segments, joints, and actuators typical of most limbs. Mechanical modeling becomes difficult and less reliable when the discrete entities interact through distributed connective tissues. In the hand, this arises when tendons from multiple muscles insert onto capsular structures around finger joints rather than directly onto individual bones and when separately controlled neuromuscular compartments of multiheaded digit muscles are loosely coupled to each other (Schieber and Santello, 2004). Distributed viscoelasticity in the skin tends to dominate the mechanics of the lightweight digits at low muscle recruitment. In the neck and trunk, modeling challenges arise when the length and/or pulling direction of one muscle depends on the position or activation of other muscles (Richmond et al., 2001). In the shoulder and hip, this occurs when muscles wrap over and around other muscles (van der Helm, 1994). In theory, these complexities can be approximated by decomposition into multiple, discrete and classical entities. In practice, most of the parameters that need to be specified to define such finite element models are unavailable.

Models of the rapidly conducting proprioceptors (myelinated nerve fiber groups I and II) are well-developed but they must be driven by perhaps uncertain musculoskeletal dynamics (muscle fascicle length and velocity for spindles and active muscle force for Golgi tendon organs). Several other types of mechanoreceptors have been described in muscles, ligaments, and joint capsules (Grigg and Hoffman, 1984), but mathematical models are not available and their functional roles are uncertain.

## Non-stationarity of Muscle Physiology

Models of the contractile properties of mammalian skeletal muscle are perhaps the most developed of all components. They were developed initially to explore reductionist models of contractile mechanisms under highly limited and generally unphysiological operating conditions (e.g., isometric or isotonic twitch or tetanus). They were extended to account fairly well for the full range of kinematic conditions and contractile properties of mammalian slow-twitch and fast-twitch muscle fibers at various physiological levels of recruitment. More recently they were extended to account for energy consumption (Tsianos et al., 2012), which may be an important cost-function used by the brain to improve performance during motor learning. What is missing are models of how the properties of muscles change over time as a result of patterns of use or disuse. The control strategies learned by the brain must anticipate or at least cope well with these changes. Understanding such changes is important when modeling is used to account for performance in pathological systems (see below).

In the short term, muscles are subject to fatigue—a reduction in force output for a given set of operating conditions (Enoka and Duchateau, 2008). Such reductions may arise from changes in the many cellular processes involved in muscle activation and deactivation. Models of muscle that are composed of computational elements that correspond to the energyconsuming processes (e.g., cross-bridge turn-over and calcium flux) should make it possible to account for those changes, but this remains to be modeled. The effort will require a great variety of experimental data, much but perhaps not all of which is already available in the literature. Many diseases and injuries are associated with disuse atrophy of muscles, which tends to increase greatly their susceptibility to fatigue, but the relative contributions of the underlying processes may be different from normal muscles.

Musculotendinous injuries account for the majority of emergency medical treatment and long-term disabilities, so there is a great deal of interest in how they occur, how they affect performance and how they heal. Structural engineering has benefitted greatly from finite-element analytical modeling (FEM) of complex structures. There have been a few attempts to develop such models for muscle, which is essentially a composite material consisting of contractile elements (muscle fibers) embedded in a structurally critical matrix (endomysial collagen) (Trotter, 1993).

Muscle is particularly responsive to exercise, which results in rapid and profound changes in key properties such as forceoutput, speed of contraction and relaxation, and resistance to fatigue. Quantitative models of these changes are becoming more feasible as the cellular mechanisms responsible for signaling and managing plasticity are starting to be revealed. Data from formal physiological experimentation and the effects of athletic training tend to emphasize covariances in physiological properties; e.g., intense but brief exercise tends to develop fast-twitch muscle fibers with larger force output, faster rise and fall dynamics, higher rates of maximal shortening (Vmax) and lower fatigue resistance than slow-twitch muscle. Mathematical models that represent explicitly the mechanisms underlying these properties are a better starting point for models of plasticity because the properties do not necessarily covary. Little systematic data are available about the rates at which the individual properties other than maximal isometric force change during normal training. For example, muscles subject to disuse atrophy produce lower force (typical of slow twitch fibers) but also have lower fatigue resistance (typical of fast twitch fibers).

The anatomy and physiology of mature neuromusculoskeletal systems are the result of myriad mechanisms that orchestrate their development (Crawford and Horowits, 2011). Most musculoskeletal models assume that the pinnation and sarcomere lengths of muscle fascicles tend to be optimized for the range of lengths that the muscles experience during normal function and that the connective tissues that support them are matched to the stresses that the contractile elements apply during those functions. The trophic factors that drive the deposition and structural properties of collagen in tendons, aponeurosis, and endomysium are starting to be understood qualitatively but quantitative models of their dynamics are not yet available. Assumptions of covariance and optimality are likely to breakdown if models are applied to pathological conditions. For example, passive tensile properties of mammalian muscle are dominated by the endomysial connective tissue, so do not necessarily covary with active tension or force-length properties that are determined by the myofilaments (Brown et al., 1996).

#### Tactile Mechanics and Sensory Transduction

Musculoskeletal systems interact with the world via contact regions composed of skin and related epidermal structures (nails, claws, teeth, hair, or fur). The mechanical compliance of skin defines what happens during object manipulation, which is essentially a series of collisions between the object and various parts of the hand; the resulting deformations of the skin define what information will be available to the CNS from tactile mechanoreceptors, which are all essentially strain gauges. The mechanical properties of these interfacial regions are starting to be captured by finite element modeling (Kumar et al., 2015). Such models can be added to classical free-body models to describe accurately the mechanical events that occur during collisions and manipulation between, for example, fingertips and objects to be grasped, but the computational load is often daunting. Modeling the tactile sensory signals that will result from myriad, independent receptors in the soft tissues remains challenging. Such sensory signals are known to be essential for dexterity, which is severely impaired when skin is anesthetized while proprioception and motor commands are left intact. The histological structure and physiological responsiveness of the receptor modalities are fairly well-understood, so it should be possible, albeit computationally challenging, to integrate model populations of cutaneous mechanoreceptors into finite element models (Kumar et al., 2015).

# Gaps in Spinal Cord Models

# Fusimotor Control of Muscle Spindles

Muscle spindles have long figured prominently in theories of sensorimotor control. The size and speed of their sensory axons, their numbers, and distributions, and their highly evolved structure and fusimotor control mechanisms all suggest that they are functionally critical. All of these elements are well-described in existing computational models (Mileusnic et al., 2006). The problem is that only sparse information is available about how their sensitivity is controlled by the fusimotor system during natural behaviors (Loeb, 1984; Prochazka, 1999; Taylor et al., 2000).

The fusimotor apparatus has undergone a huge elaboration in mammals, including specializations within various parts of the musculoskeletal system (Richmond et al., 1986). Servocontrol models of sensorimotor control originally emphasized the role of spindle afferents in the clinically prominent monosynaptic stretch reflex, but this represents only a tiny fraction of their central projections. Spindle afferents contribute to a multitude of oligosynaptic circuits in spinal and brainstem interneurons where they are combined with descending command signals. They are also the dominant source of information about posture and kinesthesia (Scott and Loeb, 1994; Gandevia, 1996), which are necessary for high-level planning and evaluation of motor strategies in cerebral cortex and cerebellum. Simplistic rules for fusimotor control have been hypothesized [e.g., alpha-gamma coactivation (Vallbo, 1974), optimal transducer programming (Loeb and Marks, 1985)] but these are speculations rather than known facts.

### Connectivity of Spinal Interneurons

Signals from the brain and sensory afferents are mixed in spinal interneurons that project to the alpha motoneurons (Pierrot-Deseilligny and Burke, 2005). As a result, the influence of brain activity on muscle contractions and movement depends substantially on the connectivity of spinal interneurons. In the past century, many neural pathways from sensory receptors (cutaneous and proprioceptive) to alpha motoneurons have been identified. These pathways include the monosynaptic Ia excitation of alpha motoneurons and polysynaptic pathways involving propriospinal, Renshaw, Ia, and Ib interneurons. The pathways were identified mainly through electrophysiology techniques by perturbing sensory signals from one muscle and investigating the timing of enhancement or depression of alpha motoneuron activity of the same muscle as well as muscles that were functional antagonists, synergists, or both (for a recent review, see Loeb, 2014). Neural pathways involving cutaneous receptors are poorly defined relative to those involving proprioceptors because it is difficult to stimulate a large group of afferents from a homogeneous type of cutaneous receptor. In electrophysiology experiments, the connectivity of cutaneous pathways is investigated typically by stimulating cutaneous nerves or patches of skin, which contain neurons from many different types of cutaneous receptors (Pierrot-Deseilligny and Burke, 2005). Through this method, it is not possible to distinguish the contribution of each type of receptor on the resulting alpha motoneuron activity. Information from these experiments, therefore, cannot be used to accurately predict the effects of physiological activity of the various cutaneous receptor afferents on interneuron or alpha motoneuron activity in the spinal cord.

Neural pathways from cutaneous receptors and proprioceptors to alpha motoneurons involving three synapses or more are not well-defined because current experimental techniques are limited. In electrophysiology experiments, if one of the neurons in the polysynaptic chain is hyperpolarized substantially via descending control, the whole chain will be invisible to the experimenter. Even if the effect on alpha motoneuron activity is observable, it will likely be highly variable and even differ in sign across subjects and experiments because descending control of each one of the neurons in the chain will likely depend strongly on the conditions of the experiment and physiological state of the subject. Furthermore, sensory afferent collaterals may activate several pathways in parallel with similar latencies, making them indistinguishable in the alpha motoneuron signal.

# Role of Presynaptic Modulation

Neurotransmitter release by the presynaptic neuron depends not only action potential rate, but also on the level of presynaptic inhibition/facilitation induced by interneurons forming axoaxonic connections with the presynaptic neurons. The activity of these interneurons are controlled by descending projections as well as by sensory afferents, but this connectivity is poorly understood (Rudomin and Schmidt, 1999). In theory, presynaptic control could allow selective gating of signals from different sensory modalities, muscles, and interneurons depending on the task. It is not clear, however, what portion of the presynaptic terminals the brain can control independently (Sirois et al., 2013). Furthermore, it is not clear how presynaptic input varies across tasks. There may be a set of presynaptic inputs that are useful for a wide range of tasks (see Fink et al., 2014), thereby reducing the number of control parameters the brain would have to learn to perform new movements.

# Status of Brain-level Models

Models of different parts of the brain tend to be abstract because of the lack of specific knowledge of the neural circuits that process signals from their input sources and the lack of knowledge of the specific neural circuits that their output projections influence. Spinal interneurons receive input from many parts of the brain, including cerebral cortex, brain stem, and tectum. The distribution of these inputs among spinal interneurons as well as the specific source locations are poorly understood. New tools for tracing and modulating connections based on genetic engineering of specific cell-types are just starting to reveal the functional relationships between brain and spinal cord circuits (Akay et al., 2014; Azim et al., 2014; Esposito et al., 2014).

# Cerebral Cortex

Voluntary motor control is usually assumed to reside in a small set of frontal lobe cortical areas. Damage to these areas in humans such as from stroke results in profound losses of such behaviors. Top-down models often attribute to cortex most or all of the control function involved in learning and executing these behaviors. Cortex seems to be necessary for learning new motor skills but may not be sufficient (see below) or even necessary for their execution (Kawai et al., 2015).

It has long been known that the cerebral cortex is organized in layers that are associated with specific types of neurons, input sources and output destinations. These could provide the substrate for bottom-up modeling. However, little is known about the local connections among neurons from different layers and even less about the connectivity among the thalamus and distant cortical columns that span the many cortical areas involved in sensorimotor control (Hooks et al., 2013; Kaneko, 2013). For this reason, circuit models tend to rely heavily on correlations of activity between a very small subset of cortical neurons either in behaving animals and/or in response to electrical stimulation (Chadderdon et al., 2014). The apparent circuits derived from these studies apply only to the small number of neurons observed and may deviate substantially from true circuitry because there are many pathways with several interneurons between recording and stimulation sites; these interneurons may block activity through some pathways. Interneuron activity depends on experimental conditions so it is likely that the net excitatory/inhibitory influence observed during the highly constrained experiments does not apply to many sensorimotor behaviors. The specific output pathways to cerebellum and subcortical structures such as the brain stem nuclei, propriospinal interneurons, and segmental interneurons are also poorly understood.

The cerebral cortex has long been an attractive candidate for computational models because of its obvious importance for learning new tasks. The rise of digital logic in the 1940s and 50s coincided with the development of electrophysiological methods to study the activity of individual neurons, leading to the compelling notion of neurons as logical and gates (McCulloch and Pitts, 1943) and learning as rules for changing the weighting of the inputs (Hebb, 1949). Models of cortical function have mostly been elaborations of this basic scheme (Marr, 1970). The problem is that the principal output cells of the cortex appear to be vastly more complex in their computational functions. Each of the numerous tiny spines that extend from their dendrites appears to function as a sophisticated temporospatial signal processor whose output gain can be individually adjusted before contributing to the all-or-none output of the neuron as a whole (Polsky et al., 2004; Jadi et al., 2012). Ambitious attempts are underway to develop computational models of cerebral cortex based on exhaustive analysis of neural connectivity (Markram, 2006), but the computational algorithm for the individual cortical neurons is now less clear than was originally assumed.

## Basal Ganglia—Thalamus

The thalamus relays and processes information from the cerebral cortex, cerebellum, brain stem, and spinal cord, and its output neurons project to various areas of the cerebral cortex. Despite its important role, the interneuron connectivity of the thalamus and associated circuitry in the basal ganglia (Bosch-Bouju et al., 2013) is poorly understood and it is commonly omitted in models of sensorimotor control. Current models of this subsystem (Nelson and Kreitzer, 2014) are still based on simplistic analogies between neurons and electronic logic gates (McCulloch and Pitts, 1943), despite data demonstrating that their circuits and even individual neurons generate complex patterns of spontaneous activity (Llinas, 1988; Nakamura et al., 2014). Computational models of individual neurons with such properties are in their infancy. Incorporating them into large-scale systems requires many more parameters and assumptions about their individual properties and connectivity patterns, as well as much greater computing power.

#### Midbrain Tectum

The tectum (superior and inferior colliculus in mammals) may play a much larger role in sensorimotor interaction with external objects than is generally acknowledged. The superior colliculus has been intensively studied and modeled for its role in directing saccadic eye movements to visual stimuli (Fecteau and Munoz, 2006). The inferior colliculus has been studied mostly as the "relay" in which sound localization information from both ears is conveyed to the auditory cortex (Slee and Young, 2014). But in fish and amphibia, the analogous subsystem controls most of the purposeful motor behaviors of the organisms, which must be fast, accurate and well-coordinated. In all vertebrates, all exteroceptive senses capable of providing localization information about an external object (vision, hearing, and touch) converge on the tectum. The tectal outputs project to brainstem nuclei and spinal cord layers that control all of the muscles required to acquire such targets, whether by saccadic gaze movements of eye and head (and auricular pinnae in most species other than primates) or reaching movements of the limbs (and jaws and tongue in many species) (Saitoh et al., 2007; Kozlov et al., 2014; Philipp and Hoffmann, 2014). The direct projections from retina to superior colliculus are known to be the source of express saccades (Munoz and Wurtz, 1992), accurately directed gaze movements that can acquire targets about twice as fast as the transcortical loop (lateral geniculate to occipital cortex and frontal eye fields) that actually projects back to the superior colliculus, which appears to control the metrics of the saccade. There have been several reports of extremely short latency corrections of reaching movements to visual targets that shift position (Gribble et al., 2002; Wilmut et al., 2006; Perfiliev et al., 2010) and in avoidance of obstacles to locomotion (Weerdesteyn et al., 2004), which would be consistent with the tectum performing the same function in limb movements as in gaze movements. If true, this would affect profoundly the assumptions made about the computational tasks of the motor and parietal cortical regions associated with reaching in extrapersonal space, perhaps the most common task for which top-down models have been constructed.

# Cerebellum and Brainstem

The cerebellum makes an excellent poster child for the problem of modeling any one of the many subsystems that subserve sensorimotor function. Its cytoarchitecture and neurophysiology is highly distinctive, relatively simple and homogeneous and better documented than any other part of the central nervous system. It is important enough to have been endowed with over half of the neurons in almost any vertebrate. Lesions to the adult cerebellum produce profound and distinctive motor deficits, yet humans born with substantial cerebellar atresia function surprisingly well (Walker, 1944). The cerebellum has attracted several detailed computational models of its function (Albus, 1975; Marr and Thach, 1991; Kawato and Gomi, 1993; Kawato and Samejima, 2007), yet it can be plausibly argued to subserve either motor (Thach, 2014) or sensory (Bower, 1997) function, to the virtual exclusion of the other.

Although the neurons of the cerebellum and their connectivity are understood better than most structures in the brain, the precise input sources and their target neurons as well as the precise destinations of the output projections are not wellknown. For example, the climbing fibers from the inferior olivary nucleus in the medulla provide a substantial portion of the input to the cerebellum, but little is known about the inputs to the inferior olive and how they are processed. Climbing fiber activity is known to be influenced by the reticular formation and red nucleus of the brain stem, several areas from the cerebral cortex, as well as by signals from the periphery related to proprioception and touch (Brown et al., 1977; Stecina et al., 2013). Specific sources of these signals, how they are processed in the inferior olive, and how the processed output affects cerebellar function is not known. Another major source of input to the cerebellum is the mossy fiber projections from the lateral reticular nucleus of the reticular formation. There is a large degree of convergence of afferent activity ascending from the spinal cord with activity descending from the cerebral cortex, tectum, red nucleus and other parts of the brain stem (Alstermark and Ekerot, 2013). The precise source of these inputs as well as the interneuronal circuits that integrate them is largely unknown.

The brain stem possesses many subdivisions including the inferior olive, vestibular nuclei, reticular formation and red nucleus with distinct set of inputs, outputs, and interneurons. Although they have an important functional role for even simple reach and grasp movements (Alstermark and Isa, 2012), the connectivity of their neural circuits as well as the precise location of their inputs and outputs are poorly understood (Kennedy, 1990; Kuchler et al., 2002). Much of the output of cerebellum is targeted to these structures, which also receive the proprioceptive and vestibular sensory signals required to coordinate movements across the entire body. Most studies of motor behavior focus on the "prime mover" muscles most closely associated with the task, but the freely mobile, multiarticulated structures of most terrestrial animals are actually quite difficult to stabilize because of intersegmental Coriolis forces (Lackner and DiZio, 2000). For example, human subjects asked to flex and extend their elbow tend to keep the wrist still without thinking, whereas a similar weight connected by a hinge to a stick being wagged up and down would flop uncontrollably. This problem could be solved by simply cocontracting the muscles of the wrist to stiffen the joint, but that would be energetically inefficient. Instead, subjects deftly time the activation of the many different wrist muscles to cancel the rapidly changing Coriolis forces. This problem is greatly magnified when trying to maintain balance of the whole body during any rapid movement of any body part, a problem that is constantly changing as the musculoskeletal system develops and ages. Perhaps the best studied model for adaptive tuning of stabilizing reflex gains is the vestibulo-ocular reflex, which involves changes in brainstem nuclei mediated by cerebellar plasticity (Kawato and Gomi, 1992; Raymond et al., 1996; Clopath et al., 2014). The details of these circuits and their learning rules remain contentious and it is uncertain if they generalize to other sensorimotor behaviors.

# Use of Systems Models to Demonstrate Competency

Many anatomical structures of the sensorimotor system have been investigated to characterize their inputs, outputs, and intermediate processing (although as explained above, many important knowledge gaps remain). Their complexity often makes it impossible to intuit their input-output transformation. Neural networks within these structures, for example, have a large number of neurons with various non-linear properties and substantial convergent, divergent, and recurrent connections that make it difficult to infer how they process a set of incoming signals. Computational models are therefore essential for predicting these non-intuitive interactions and provide insight into the possible transformations that a neural network can apply to a range of input signal patterns. Intuition is even less useful for predicting the collective input-output transformation of a system of interconnected neural networks formed by distinct anatomical structures in the nervous system such as those shown in **Figure 1**. Understanding such complex interactions is necessary for determining the relative roles of these anatomical structures in sensorimotor control and ultimately unraveling the mechanisms involved.

The sparse knowledge of the anatomy and physiology underlying the sensorimotor system precludes a complete understanding of the mechanism of sensorimotor control. This knowledge, however, can be exploited to understand the competency of the known properties and suggest experimental investigations for furthering our understanding and updating the models. It can be tested, for example, if a set of known properties of the sensorimotor system is sufficient for generating specific behaviors. If the properties are sufficient for generating a particular behavior then this would suggest that these properties may play a significant role. It has been shown recently that the known spinal circuits described above are sufficient for generating the muscle dynamics of wrist (Raphael et al., 2010) and arm movement (Tsianos et al., 2011, 2014), which has been traditionally assumed to arise largely from commands issued by the motor cortex. These results emphasize the possibility that spinal circuits can make large contributions to voluntary movement and encourage further experimental testing. It has also been shown that descending commands to models of known spinal circuitry can be linearly interpolated to generate intermediate movements (Tsianos et al., 2014). Such simple interpolation along with the modeled musculoskeletal system and spinal circuitry properties were sufficient to reproduce the extent of learning generalization observed experimentally for similar tasks. This result suggests that much of the generalization of learning observed experimentally may arise from simple combinations of learned voluntary commands rather than through the use of internal inverse models of the sensorimotor system. This is consistent with the tendency of human subjects to adapt to changes in the musculoskeletal system by relatively simple scaling of the original motor programs rather than computation of new programs that offered superior performance (de Rugy et al., 2012a,b).

# Use of Systems Models to Prove Insufficiency

If modeled properties of the sensorimotor system are not sufficient to explain a particular behavior, then these properties must interact with other biological properties that are either not modeled or not known. Other known properties can be added to the model incrementally to test which configurations are sufficient for reproducing the desired behavior. For example, activation of muscles and feedback through spinal circuits contribute to stable movements due to the relatively short time delays involved; however, the contribution of spinal circuits is limited for whole body tasks that require muscle coordination between distant body parts. Spinal circuits coordinate activity among a limited number of adjacent muscles and joints because inputs to a given spinal interneuron originate from only a few adjacent spinal segments. There are many relatively fast circuits involving anatomical structures rostral to the spinal cord such as the brain stem (Esposito et al., 2014), tectum (Philipp and Hoffmann, 2014), cerebellum (Azim et al., 2014) and even sensorimotor cortex (Scott, 2004, 2012; Nashed et al., 2015) that receive afferent signals from throughout the body and may therefore contribute to the dynamics and stability of such movement. Known circuits from these anatomical structures can be incorporated in models to test their sufficiency, which would depend on their specific connectivity as well as the range of input sources and their extent of convergence. If multiple anatomical structures appear sufficient, then the modeled task could be varied systematically to search for situations where each one becomes uniquely competent. This would indicate the types of sensorimotor tasks and conditions that each anatomical structure might subserve in the biological system. If known properties of a given anatomical structure cannot account for the behavior, it could be that the model for the

structure is inadequate or that other structures are required. Such negative results are particularly valuable in identifying opportunities to advance understanding of the system as a whole.

# General Obstacles to Successful Modeling

It is always easy to end a scientific paper with a call for more experiments to provide more data to fill in the missing pieces of knowledge. This brief review has pointed out many places where basic knowledge about neural connectivity is insufficient to permit reductionist modeling. But there are also many places where the connectivity is relatively well-known and detailed models have been constructed, only to expose contentious disagreements about what role the structure actually plays in a given behavior.

As David Marr pointed out, models must start with a top-level theory of computation—a division of function into a sequence or hierarchy of tasks (Marr, 1982). Only then can one contemplate the computational algorithms that might subserve each of those tasks and, below that, the machinery that performs each algorithm. The physiology and connectivity of individual neurons provides information about the machinery level, leaving the modeler to guess about the tasks and the algorithms.

It used to be feasible for a researcher to become familiar with most of the literature about most of the CNS and to construct a systems level model of behaviors that considered both its phylogenetic and ontological origins; for example, (Ayres, 1975). Most researchers now spend a lifetime learning the details of and formulating hypotheses about one or two of the many subsystems that somehow contribute to sensorimotor behaviors. They tend naturally to assume that the subsystem that they are studying is primarily responsible for the behaviors that they employ in their experimental designs. Their models of "their" subsystem must nevertheless be reconciled with some current dogma about the role of the other subsystems, as promulgated by other researchers with equally narrow and parochial views.

"It was six men of Indostan To learning much inclined, Who went to see the Elephant (Though all of them were blind), That each by observation Might satisfy his mind...

And so these men of Indostan Disputed loud and long, Each in his own opinion Exceeding stiff and strong, Though each was partly in the right, And all were in the wrong!"

—From "Blind Men and the Elephant" by John Godfrey Saxe (1816–1887).

# References

for normal and parkinsonian conditions. Front. Comput. Neurosci. 7:163. doi: 10.3389/fncom.2013.00163

and three-dimensional finite element models for tactile sensation studies. J. Biomech. Eng. 137, 061002. doi: 10.1115/1.4029985

recordings from hindlimb muscle spindles. J. Physiol. 522, 515–532. doi: 10.1111/j.1469-7793.2000.t01-3-00515.x

**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.

Copyright © 2015 Loeb and Tsianos. 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.

# Direct and indirect spino-cerebellar pathways: shared ideas but different functions in motor control

Juan Jiang1† , Eiman Azim2† , Carl-Fredrik Ekerot 3† and Bror Alstermark <sup>1</sup> \* †

<sup>1</sup> Department of Integrative Medical Biology, Section of Physiology, Umeå University, Umeå, Sweden, <sup>2</sup> Departments of Neuroscience and Biochemistry and Molecular Biophysics, Howard Hughes Medical Institute, Kavli Institute for Brain Science, Mortimer B. Zuckerman Mind Brain Behavior Institute, Columbia University, New York, NY, USA, <sup>3</sup> Department of Experimental Medical Science, University of Lund, Lund, Sweden

The impressive precision of mammalian limb movements relies on internal feedback pathways that convey information about ongoing motor output to cerebellar circuits. The spino-cerebellar tracts (SCT) in the cervical, thoracic and lumbar spinal cord have long been considered canonical neural substrates for the conveyance of internal feedback signals. Here we consider the distinct features of an indirect spino-cerebellar route, via the brainstem lateral reticular nucleus (LRN), and the implications of this pre-cerebellar "detour" for the execution and evolution of limb motor control. Both direct and indirect spino-cerebellar pathways signal spinal interneuronal activity to the cerebellum during movements, but evidence suggests that direct SCT neurons are mainly modulated by rhythmic activity, whereas the LRN also receives information from systems active during postural adjustment, reaching and grasping. Thus, while direct and indirect spinocerebellar circuits can both be regarded as internal copy pathways, it seems likely that the direct system is principally dedicated to rhythmic motor acts like locomotion, while the indirect system also provides a means of pre-cerebellar integration relevant to the execution and coordination of dexterous limb movements.

#### Edited by:

Ning Lan, Shanghai Jiao Tong University, China

#### Reviewed by:

Irina N. Beloozerova, Barrow Neurological Institute, USA Vincent C. K. Cheung, The Chinese University of Hong Kong, Hong Kong

#### \*Correspondence:

Bror Alstermark, Department of Integrative Medical Biology, Section of Physiology, Umeå University, Linneus väg 9, 901 87 Umeå, Sweden<EMAIL_ADDRESS>
†These authors have contributed equally to this work.

> Received: 13 March 2015 Accepted: 01 June 2015 Published: 06 July 2015

#### Citation:

Jiang J, Azim E, Ekerot C-F and Alstermark B (2015) Direct and indirect spino-cerebellar pathways: shared ideas but different functions in motor control. Front. Comput. Neurosci. 9:75. doi: 10.3389/fncom.2015.00075 Keywords: lateral reticular nucleus (LRN), spino-cerebellar pathways, spino-LRN-cerebellar pathways, internal feedback, motor control

# Introduction

Cerebellar circuits are of major importance in the control of movements, providing a neural basis for pattern recognition and motor behavioral correction and adaptation (Ito, 2006). These contributions to motor control depend on specific mossy fiber and climbing fiber cerebellar inputs that convey information about both ongoing motor output and external sensory events (Ito, 1984; Dean et al., 2010). In this paper, we focus on the organization of mossy fiber systems and, more specifically, we delineate two classes of spino-cerebellar pathways: direct spino-cerebellar projections, and indirect pathways via the lateral reticular nucleus (LRN; referred to as the spino-LRN-cerebellar pathway), as illustrated schematically in **Figure 1**.

Spino-cerebellar pathways have been implicated in the transmission of information about external events from various sensory modalities (cf. review, Stecina et al., 2013), the cancellation of reafferent sensory signals during self-generated movements (Hantman and Jessell, 2010), and the conveyance of internal copies of motor commands for rapid motor prediction and correction (Lundberg, 1971; Arshavsky et al., 1972, 1978; Alstermark and Isa, 2012;

Fedirchuk et al., 2013; Azim and Alstermark, 2015). However, little is known about the organizational and functional logic underlying the existence of two separate systems for conveying spinal signals to the cerebellum. By comparing the phylogeny, anatomy, genetic identities and functional organization of direct and indirect spino-cerebellar circuits, we highlight key similarities and differences between these pathways and discuss principal questions that remain.

# Phylogeny

A phylogenetic comparison of cerebellar circuits has been extensively reviewed by Ito (1984). There is evidence that direct spino-cerebellar and indirect spino-LRN-cerebellar tracts coexist in teleost fish (Szabo et al., 1990; Finger, 2000), suggesting an early evolutionary divergence of these pathways. In mammals, several direct spino-cerebellar tracts (SCT) have been identified anatomically and electrophysiologically. Two of the most studied are the dorsal (DSCT) and ventral (VSCT) spino-cerebellar tracts that originate in the thoracic and lumbar spinal cord (Jankowska et al., 2011; cf. review: Stecina et al., 2013). The corresponding direct SCT for forelimb regions are the cuneocerebellar tract (CCT; Jansen and Brodal, 1954; Ekerot and Larson, 1972) in the brainstem and the rostral spino-cerebellar tract (RSCT; Oscarsson, 1965; Hirai et al., 1976) in cervical segments, respectively (**Figure 1**).

Indirect spino-LRN-cerebellar pathways have mainly been studied in the cat (Clendenin et al., 1974a,b,c,d, 1975; Matsushita and Ikeda, 1976; Ekerot, 1990a,b,c), but comparative anatomy (Walberg, 1952) has revealed the existence of the LRN in a large number of mammals including Erinaceomorpha (hedgehog), Chiroptera (bat), Rodentia (squirrel, mouse and rat), Lagomorpha (hare), Carnivora (cat, dog and seal), Cetartiodactyla (harbor porpoise), Artiodactyla (pig, cow and roe deer) and Primates (rhesus macaque and human).

Interestingly, in teleosts there are abundant axon collaterals from direct spino-cerebellar pathways to the LRN (Szabo et al., 1990), whereas in cats the DSCT does not provide collateral excitation to LRN neurons (Ekerot and Oscarsson, 1975). These phylogenetic differences suggest that direct spino-cerebellar and indirect spino-LRN-cerebellar pathways may have originated as cooperative systems, which became progressively separated as more advanced motor repertories evolved.

# Anatomy

As shown in **Figure 1**, direct spino-cerebellar and indirect spino-LRN-cerebellar pathways originate in cervical, thoracic and lumbar spinal segments, as well as in the brainstem (for review, cf. Alstermark and Ekerot, 2013; Pivetta et al., 2014). Within the direct and indirect classes, subpopulations with ipsilateral, contralateral and bilateral projections have been identified (cf. reviews: Alstermark and Ekerot, 2013; Stecina et al., 2013). The ultimate mossy fiber terminations of these pathways in the cerebellar cortex are found mainly in the vermal and paravermal regions of the anterior and posterior lobes, as well as in the paramedian lobe. The location of ascending axonal projections in the white matter of the spinal cord and the pattern of mossy fiber termination zones within the cerebellar cortex differ across individual systems, but broad comparison of direct and indirect pathways to each other reveals no clear differences (cf. review Ito, 1984). Thus, at least at the gross anatomical level, direct spino-cerebellar and indirect spino-LRN-cerebellar pathways target overlapping cerebellar cortical circuits.

# Genetic Identities

The genetic delineation of neuronal subtypes has complemented classical anatomical and electrophysiological characterization of spinal circuits, and has provided a means for selective manipulation and functional dissection of these pathways (Goulding, 2009). While the molecular identities of each of the direct spino-cerebellar systems are yet to be fully defined, studies in mice have revealed that a population of dorsally-derived spinal interneurons that express the transcription factor Math1 give rise to multiple spino-cerebellar pathways (Bermingham et al., 2001). Moreover, DSCT neurons in Clarke's column have been shown to selectively express the neurotrophic factor Gdnf (Hantman and Jessell, 2010).

Indirect spino-LRN-cerebellar pathways, and the cervical propriospinal neuron (PN) system in particular, have been the subject of much recent genetic scrutiny. A prominent population of excitatory PNs involved in goal-directed reaching movements was identified within the Chx10-expressing V2a interneuron class (Azim et al., 2014); notably, only cervical but not lumbar V2a interneurons project to the LRN, indicating that indirect LRN-cerebellar pathways originating in the lumbar cord have distinct genetic identities. In zebrafish, a subset of V2a spinal interneurons send ascending projections to the hindbrain (Menelaou et al., 2014), suggesting that the V2a interneuron class establishes an evolutionarily conserved circuit for the conveyance of motor signals to supraspinal regions. Moreover, recent genetic and viral labeling studies in mice have revealed that in addition to V2a interneurons, several classes of molecularly defined excitatory and inhibitory cervical spinal interneurons project to the LRN (Pivetta et al., 2014), suggesting that other indirect spino-cerebellar pathways can be dissected genetically along similar lines.

# Functional Organization

It has been well documented that both direct spino-cerebellar and indirect spino-LRN-cerebellar pathways convey information related to ongoing rhythmic movements, including locomotion, scratching and respiration (cf. references in reviews by Ito, 1984; Alstermark and Ekerot, 2013; Stecina et al., 2013). Moreover, it has been proposed that the VSCT (Lundberg and Weight, 1971) and DSCT (Hantman and Jessell, 2010) monitor the excitability of spinal interneurons. Interestingly, whereas the VSCT (Fedirchuk et al., 2013) and DSCT (Stecina et al., 2013) signal mainly during the flexion phase, spino-LRN-cerebellar pathways are active throughout the entire cycle of flexion and extension (cf review by Alstermark and Ekerot, 2013), suggesting that indirect pathways convey a broader range of motor signals.

Another major difference in the functional organization of direct and indirect cerebellar pathways is that the four subsystems in the indirect spino-LRN-cerebellar pathway originating in the cervical spinal cord and brainstem (**Figure 1**) may be dedicated to more than just rhythmic movements (Alstermark and Ekerot, 2013). These subsystems, by monitoring the excitability of spinal interneurons, could signal information about posture (bilateral ventral flexor reflex tract; bVFRT), reaching (C3-C4 propriospinal system; PN), grasping (ipsilateral forelimb tract; iFT) and jaw opening (dorsal funiculus-trigeminal tract; DF-Trig), and their convergence in the LRN might enable the coordination of these separate motor actions into coherent and smooth movements (Alstermark and Ekerot, 2013, 2015).

Among these indirect systems, the function of C3-C4 PNs has been investigated extensively in the cat, monkey, human and recently in the mouse (Alstermark and Isa, 2012; Azim et al., 2014). These studies have shown that PNs mediate motor commands for reaching by directly modulating the activity of forelimb-innervating motor neurons, while also conveying copies of these motor commands, via axon collaterals, to the LRN. Genetic manipulation of PNs in the mouse has revealed that this internal copy pathway recruits a cerebellarmotor feedback loop, providing a plausible neural substrate for the rapid updating and correction of ongoing forelimb motor output (Azim et al., 2014; Azim and Alstermark, 2015). The current lack of selective genetic access to other spino-LRN-cerebellar pathways has precluded similar exploration of their behavioral functions, yet evidence suggests that the cervical bVFRT, iFT and PN systems provide both discrete and convergent internal feedback signals to LRN-cerebellar circuits (Alstermark and Ekerot, 2013; Pivetta et al., 2014; Huma and Maxwell, 2015), potentially enabling the coordination of forelimb and postural motor control. A companion article discusses the extensive convergence of projections from these distinct systems in the LRN, providing a pre-cerebellar center for the integration of spinal signals and their modulation by descending motor cortical pathways (Alstermark and Ekerot, 2015).

# Open Questions and Future Directions

1. How do descending motor pathways modulate the direct and indirect spino-cerebellar tracts? Thus far, only the descending inputs to the cervical PN system have been investigated systematically (Alstermark and Lundberg, 1992; Alstermark and Isa, 2012; Azim et al., 2014). The convergence of descending pathways onto cervical PNs suggests a role for these neurons in integrating motor command signals and conveying copies of this information to LRN-cerebellar circuits. A better understanding of the descending inputs onto other direct and indirect cerebellar pathways should help to clarify their potential contributions to voluntary movements.

# References

particular suggests that these systems may have evolved in concert with the increasing complexity of dexterous forelimb movements. The identification of unique genetic markers for each of these pathways should offer a means to access and manipulate these circuits selectively, providing the experimental resolution needed to characterize their discrete contributions to motor behavior (Azim et al., 2014; Azim and Alstermark, 2015).

6. There is growing interest in applying computational neurobiology approaches to understanding the molecular and genetic mechanisms that may contribute to spinocerebellar ataxia (cf. review by Brown et al., 2015), and models devoted to the role of internal feedback more generally have explored various neural circuits in the cortex, brainstem and spinal cord (cf. review by Azim and Alstermark, 2015). Regarding spinocerebellar pathways, a hypothesis has been forwarded on their role in the multi-dimensional integration of sensorimotor information (Spanne and Jörntell, 2013). However, in this model, direct spino-cerebellar and indirect spino-LRN-cerebellar pathways are grouped together. A new hypothesis has recently been proposed that focuses specifically on the role of indirect spino-LRN-cerebellar pathways (Alstermark and Ekerot, 2013). Future modeling approaches, informed by the experimental work described above, should provide greater insight into the discrete functions of direct and indirect spino-cerebellar systems.

# Acknowledgments

This work was supported by grants from the Swedish Research Council and Umeå University to BA and by a National Institutes of Health (NIH) K99 award (NS088193) to EA.

**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.

Copyright © 2015 Jiang, Azim, Ekerot and Alstermark. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution and 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.

# Subject-specific computational modeling of DBS in the PPTg area

Laura M. Zitella<sup>1</sup> , Benjamin A. Teplitzky <sup>1</sup> , Paul Yager <sup>2</sup> , Heather M. Hudson<sup>2</sup> , Katelynn Brintz <sup>1</sup> , Yuval Duchin<sup>3</sup> , Noam Harel <sup>3</sup> , Jerrold L. Vitek <sup>2</sup> , Kenneth B. Baker <sup>2</sup> and Matthew D. Johnson1, 4 \*

*<sup>1</sup> Department of Biomedical Engineering, University of Minnesota, Minneapolis, MN, USA, <sup>2</sup> Department of Neurology, University of Minnesota, Minneapolis, MN, USA, <sup>3</sup> Center for Magnetic Resonance Research, University of Minnesota, Minneapolis, MN, USA, <sup>4</sup> Institute for Translational Neuroscience, University of Minnesota, Minneapolis, MN, USA*

Deep brain stimulation (DBS) in the pedunculopontine tegmental nucleus (PPTg) has been proposed to alleviate medically intractable gait difficulties associated with Parkinson's disease. Clinical trials have shown somewhat variable outcomes, stemming in part from surgical targeting variability, modulating fiber pathways implicated in side effects, and a general lack of mechanistic understanding of DBS in this brain region. Subject-specific computational models of DBS are a promising tool to investigate the underlying therapy and side effects. In this study, a parkinsonian rhesus macaque was implanted unilaterally with an 8-contact DBS lead in the PPTg region. Fiber tracts adjacent to PPTg, including the oculomotor nerve, central tegmental tract, and superior cerebellar peduncle, were reconstructed from a combination of pre-implant 7T MRI, post-implant CT, and post-mortem histology. These structures were populated with axon models and coupled with a finite element model simulating the voltage distribution in the surrounding neural tissue during stimulation. This study introduces two empirical approaches to evaluate model parameters. First, incremental monopolar cathodic stimulation (20 Hz, 90µs pulse width) was evaluated for each electrode, during which a right eyelid flutter was observed at the proximal four contacts (−1.0 to −1.4 mA). These current amplitudes followed closely with model predicted activation of the oculomotor nerve when assuming an anisotropic conduction medium. Second, PET imaging was collected OFF-DBS and twice during DBS (two different contacts), which supported the model predicted activation of the central tegmental tract and superior cerebellar peduncle. Together, subject-specific models provide a framework to more precisely predict pathways modulated by DBS.

Keywords: pedunculopontine nucleus, deep brain stimulation, Parkinson's disease, non-human primate, finite element, diffusion tensor

# Introduction

Gait and balance difficulties in Parkinson's disease can be especially debilitating since they increase the risk of falling. For some patients, these symptoms are resistant or poorly managed by levodopa treatment and typical targets of deep brain stimulation (DBS) including the subthalamic nucleus (STN) and internal segment of the globus pallidus. In contrast, low-frequency electrical stimulation delivered within or near the pedunculopontine tegmental nucleus (PPTg), a component of the

#### Edited by:

*Ning Lan, Shanghai Jiao Tong University, China*

#### Reviewed by:

*Frank Rattay, Vienna University of Technology, Austria Ursula Van Rienen, University of Rostock, Germany*

#### \*Correspondence:

*Matthew D. Johnson, Department of Biomedical Engineering, University of Minnesota, 7-105 Nils Hasselmo Hall, 312 Church Street SE, Minneapolis, MN 55455, USA<EMAIL_ADDRESS>
> Received: *13 March 2015* Accepted: *02 July 2015* Published: *14 July 2015*

#### Citation:

*Zitella LM, Teplitzky BA, Yager P, Hudson HM, Brintz K, Duchin Y, Harel N, Vitek JL, Baker KB and Johnson MD (2015) Subject-specific computational modeling of DBS in the PPTg area. Front. Comput. Neurosci. 9:93. doi: 10.3389/fncom.2015.00093* mesencephalic locomotor region (MLR) of the brainstem, has shown promising results (Plaha and Gill, 2005; Schrader et al., 2013; Mazzone et al., 2014). Clinical outcomes, however, have varied from patient to patient across these studies, due in part to variation in surgical targeting, uncertainty in the therapeutic target, and the likely modulation of highly excitable, side effect inducing fiber pathways (Nowak and Bullier, 1998) outside the MLR.

A previous computational modeling study showed that clinical outcomes of DBS within the PPTg area are likely to be highly dependent upon lead position and stimulation settings (Zitella et al., 2013). For instance, the superior cerebellar peduncle passes through the PPTg en route from the deep cerebellar nuclei to the red nucleus and cerebellar receiving area of the motor thalamus by means of a decussation just rostral to PPTg. At present, how selective activation of this pathway affects freezing of gait is unknown, though stimulation of the cerebellothalamo-cortical circuit has been postulated to be beneficial for gait (Fournier-Gosselin et al., 2013). DBS in the PPTg area may also modulate medial fiber tracts such as the medial longitudinal fasciculus (MLF) or the oculomotor nerve (ON). Side effects from activation of either of these fiber tracts would be expected to affect the eyes and eyelids. Neuronal activation volumes that extend lateral of PPTg may include the medial lemniscus (ML) and lateral lemniscus (LL) and lead to paresthesias (Murata et al., 2003) and changes in auditory perception (Lim et al., 2008), respectively. Further, spread of current rostral to the PPTg may modulate the central tegmental tract (CTG), which rises from the nucleus solitarius and carries gustatory input to the ventral posteromedial nucleus of the thalamus (VPM). There is some evidence that activation of this tract may result in palatal myoclonus (Matsuo and Ajax, 1979).

While previous subject-specific computational models of DBS have been developed for the STN (Miocinovic et al., 2006; Chaturvedi et al., 2010; Butson et al., 2011), tailoring models to the PPTg area has been limited because of the poor image contrast within brainstem with standard MR scanner technology. In recent years, however, advances in structural imaging have made visualizing fiber tracts within the brainstem more readily available. In this study, we leverage susceptibility-weighted imaging (SWI) and diffusion-weighted imaging (DWI) to create subject-specific computational models of PPTg-DBS, which can predict activation of individual fiber tracts within the brainstem for any given DBS setting. In order for these models to be informative for clinicians, the models must provide accurate predictions of neuronal activation. The challenge becomes defining behavioral or functional outcome measures to confirm or otherwise modify the selection of model parameters including tissue conductance anisotropy and inhomogeneity, cellular morphology, axonal diameter, and ion channel kinetics among others. Here, we propose two approaches in the context of PPTg-DBS, namely eliciting an oculomotor side effect and performing DBS within the context of positron emission tomography (PET) imaging.

To examine the pathways modulated in the PPTg area during PPTg-DBS, a subject-specific computational model was developed. In this study, the models were used to (1) investigate the effects of using tissue conductance anisotropy within the brainstem based on diffusion-weighted imaging, and (2) perform model parameter sweeps to determine PPTg-DBS model sensitivity. The models were evaluated with varying axon diameter, conductivity values, and DBS lead location, and then compared against behavioral and PET imaging results.

# Materials and Methods

# Subject

Two female rhesus macaque monkeys (macaca mulatta, Monkey L and Monkey P) were used in this study. All procedures were approved by the Institutional Animal Care and Use Committee of the University of Minnesota and complied with United States Public Health Service policy on the humane care and use of laboratory animals. The animals were housed individually with environmental enrichment, provided with water ad libitum, and given a range of food options including fresh fruit and vegetables. All efforts were made to provide good care and alleviate any discomfort for the animals during the study.

Pre-operative 7T MRI was acquired at the Center for Magnetic Resonance Research (CMRR) at the University of Minnesota using a passively shielded 7T magnet (Magnex Scientific) for both animals. During the imaging sessions, the animals were anesthetized with isoflurane (2.5%) and monitored for depth of anesthesia. Susceptibility-weighted imaging was acquired with a 3D flow-compensated gradient echo sequence at 0.4 mm isotropic resolution using a field of view (FOV) of 128 × 96 × 48 mm<sup>3</sup> . Diffusion-weighted images (b-value = 1500 s/mm<sup>2</sup> ) were acquired with diffusion gradients applied along 110 uniformly distributed directions using a 128 × 84 × 99 mm<sup>3</sup> FOV (1 mm isotropic resolution). The 3D tensors were calculated as ellipsoidal functions, to identify the orientation of maximum value (Barmpoutis et al., 2009; Barmpoutis and Vemuri, 2010).

In Monkey L, a cranial chamber was mounted on the head to facilitate implantation of the DBS lead, as described previously (Elder et al., 2005). The high-field imaging, along with results from electrophysiological microelectrode mapping of the PPTg area, were superimposed in Monkey Cicerone (Miocinovic et al., 2007) to define a trajectory for unilateral implantation of a scaled-down version of a human DBS lead (2F diameter, 8 annular electrode contacts: 0.5 mm height, 0.25 mm spacing) (NuMed, Hopkinton, NY) in the region of the PPTg (right hemisphere). Following lead implantation, a post-operative CT scan was performed under Ketamine and Dexdomitor anesthesia to visualize the implantation trajectory and depth in Monkey Cicerone. The preoperative SWI was coregistered with the postoperative CT to determine the DBS lead location relative to nuclei and fiber tracts within the brainstem. After instrumentation with the chamber and the DBS lead, the subject was rendered parkinsonian with systemic injections of 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP).

At the conclusion of the study, both animals were deeply anesthetized with sodium pentobarbital and perfused with a fixative solution containing 4% paraformaldehyde, consistent with the recommendations of the Panel of Euthanasia of the American Veterinary Medical Association. After fixation, the

brain was removed, blocked, and cryoprotected in 15% sucrose in phosphate buffered solution. Coronal sections (50µm) were cut using a freezing microtome and labeled for Nissl. In the case of Monkey L, the DBS lead trajectory was again reconstructed from these histological images using Mimics (Materialise, Leuven, Belgium), which confirmed the in vivo estimation of the DBS lead trajectory that had been generated from co-registration of the pre-implant SWI with the post-implant CT.

# Axonal Model Morphologies

Several imaging-based tools were used to reconstruct the threedimensional morphologies of the PPTg area for use in the subject-specific finite element model and multi-compartment neuron model simulations (**Figures 1A,B**). The SWI volume was aligned in anterior commissure to posterior commissure (AC-PC) space with Analyze (AnalyzeDirect, Overland Park, KS), and then resliced to generate images that matched atlas plates from a rhesus macaque brain atlas (Paxinos, 2009). A nonlinear affine atlas registration algorithm based on a moving least squares fit applied to each image voxel was used to identify the borders of the PPTg, CTG, and ON. The algorithm involved the selection of analogous control points placed on each MR image and each corresponding atlas plate. Contours were traced from these warped atlas reconstructions in Rhinoceros (McNeel, Seattle, WA) and lofted into 3D surfaces using a nonuniform rational B-spline modeling approach. Results from the warping algorithm were aligned to the DWI using FLIRT, which provided anatomical context to guide the placement of seed points for probabilistic diffusion tractography calculations in FSL (Jenkinson et al., 2012). To identify the SCP tract, seed points were placed in the cerebellar outflow tract caudal to the decussation, with waypoints defined at the decussation of SCP and the entire contralateral thalamus. Manual thresholding of the output of probtrackx was performed in Amira (FEI, Hillsboro, OR) to produce the final tract geometry.

## Finite Element Models (FEM)

A finite element model of the DBS lead and surrounding neural tissue was created in COMSOL (**Figure 1C**). The variable resolution tetrahedral mesh was constructed with a maximum element size of 0.2 mm for the electrodes and 1.6 mm for the lead, encapsulation layer, and neural tissue. The final mesh consisted of 447280 elements with a finer resolution near the electrode-tissue interface. A point current source was modeled at the geometric center of each electrode, and the entire lead was surrounded by a 250µm encapsulation layer with a conductivity of 0.18 S/m. A 20 mm radius sphere around the electrode represented the neural tissue. A cylinder was placed on the edge of the sphere to represent the cranial chamber, and the chamber outer surface was assigned as ground (**Figure 1C**). For the isotropic FEM, conductivity of the neural tissue was homogeneous, 0.3 S/m. For the anisotropic FEM, the neural tissue conductivity was calculated from the 6-direction DTI tensors, based on an estimated linear relationship between the conductivity (σ) and diffusion tensor eigenvalues (Tuch et al., 2001).

$$
\sigma = s \ast D \tag{1}
$$

location and grounded chamber. (D) Electric potential isosurfaces for the anisotropic and isotropic model.

where s was set to 0.844 and D represents the diffusion tensor eigenvalues. These conductivity matrices (**Figures 2A,B**) were imported into COMSOL and interpolated onto the mesh using the nearest neighbors function. Voltage distribution during monopolar stimulation through each electrode were solved for using anisotropic and isotropic conductivity models with COMSOL (**Figure 1D**) (Schmidt and Van Rienen, 2012). The extracellular voltage predicted from the FEM solution was then interpolated at each nodal compartment along each multi-compartment axon model.

# Biophysical Modeling of DBS in the PPTg Area

SCP, ON, and CTG were each randomly populated with multi-compartment cable models of 1000 myelinated axons, ranging in diameter from 2 to 8.7µm (McIntyre et al., 2004a; Miocinovic et al., 2006; Johnson and McIntyre, 2008; Birdno et al., 2011). The axonal models included nodes of Ranvier, paranodal, and intermodal segments as well as a myelin sheath (McIntyre et al., 2002). Nodal compartments were given biophysical mechanisms related to a nonlinear fast Na<sup>+</sup> channel, persistent sodium channel, slow K<sup>+</sup> channel, and a leakage current. The paranodal compartments were instantiated with a slow K<sup>+</sup> current. Both intermodal and myelin compartments had only passive mechanisms, including a membrane capacitance (Cm = 2 µF/cm<sup>2</sup> ) and axoplasmic resistivity (70 -cm). Stimulus thresholds for evoking action potentials within the modeled axons were estimated in the NEURON v7.3 programming environment (Hines and

FIGURE 2 | Comparison of conductivity and diffusion tensors between Monkey L (top) and Monkey P (bottom). (A) The calculated conductivity, σxx, is shown for select coronal slices. (B) The distribution of conductivity values calculated from the primary, secondary, and tertiary eigenvalues for the entire brain (top) and the segmented brainstem (bottom). (C) The fractional anisotropy for a select brainstem slice (left), compared to a corresponding T1 slice (right). The diffusion tensors are plotted as spherical functions and overlaid on the fractional anisotropy. The orientations (dorsal-caudal, anterior-posterior, medial-lateral) are represented as RGB color components (i.e., R, G, and B, respectively).

Carnevale, 1997). Stimulus perturbations were inserted using the NEURON mechanism, extracellular, for each axonal compartment. Stimulus pulses (90µs pulse width) were delivered at a rate of 20 Hz. Activation was defined by the lowest amplitude to elicit one or more action potentials within 1–3 ms following each stimulus pulse for at least 8 of 10 stimulus pulses.

## Motor Side-Effects of DBS in the PPTg Area

A monopolar review was conducted to determine the electrical stimulation amplitude thresholds for eliciting an overt motor side effect, which for this lead implantation was found to be a right eyelid flutter. In this case, stimulation was delivered through an externalized current-driven pulse generator (Precision, Boston Scientific) with a cathode applied individually to each contact and a return set to the cranial chamber. The stimuli were delivered as a 20 Hz train of 90µs pulses in 0.1 mA increments until the eyelid flutter was observed visually by the investigators, as shown in **Table 1**.

## PET Analysis

Approximately 1 year after implantation of the DBS lead in Monkey L, PET/CT was collected using a Siemens Biograph 64

#### TABLE 1 | Motor side-effect thresholds.

slice scanner on three different days within an 8-day period. A full 24 h before each scan, the subject was withheld from any stimulation or medication. The subject was fasted beginning at 1700 the night prior to the scan (Garraux et al., 2011), with fasting blood glucose verified the morning of the scan. Thereafter, the proximal end of the DBS lead was connected to the external pulse generator and a single 8 mCi dose of 18- FDG was administered intramuscularly. As listed in **Table 2**, immediately following this injection either 20 Hz PPTg-DBS was applied at one of the two contacts of the DBS lead or no DBS was applied (baseline scan). The subject sat still in a quiet, familiar environment without stimuli while this

TABLE 2 | DBS conditions during each PET scan.

treatment was administered. After 45 min, general anesthesia was induced using ketamine (10 mg/kg) and diazepam (0.5 mg/kg) (Oguchi et al., 1982; Wyckhuys et al., 2014), and the subject was moved to the scanner for imaging. Reconstruction yielded a voxel size of 1.018 × 1.018 × 2 mm. The subject was then released to an isolation room until radioactivity was undetectable.

The PET images were analyzed with methods similar to those described previously (Ponto et al., 2014; Wyckhuys et al., 2014). Each non-contrast CT scan was transformed by rigid, manual co-registration (Jena et al., 2014) to align with a standard MRI template, INIA19 Macaca mulatta (Rohlfing et al., 2012), with the resultant transformation values individually applied to the PET image from the corresponding data acquisition session to align them to the normalized space. Finally, a preoperative MRI, taken before the cephalic hardware and DBS lead had been placed, was aligned with the INIA19 MRI to verify fit. The 1210-MRIderived-VOIs of the INIA19 template, were used within PMOD software (PMOD Technologies, Zurich, Switzerland) to compile statistical measures. Mean voxel value, standard deviation, and number of voxels were collected as standard uptake values (SUV) (Garraux et al., 2011). All scans were scaled so that the left occipital white matter uptake was equivalent. Volumes of interest were then consolidated to yield larger brain structures, decreasing the resolution from 1210 volumes of interest into a manageable grouping for analysis. The scaled images had a two tailed t-test performed at each of the VOIs, with a noncorrected p-value (alpha = 0.05).  These brain region SUVs are reported below along with respective T scores and p-values. Due to the use of a single subject, no CT transformation based image attenuation correction was performed on the PET scan results.

# Results

Accurate prediction of therapeutic outcomes by computational neuron models of DBS targeting the MLR of the brainstem will have strong clinical value for freezing of gait in Parkinson's disease patients with these implants (Zitella et al., 2013). Two major challenges remain in rendering these computational models more predictive in power: (1) making them subject specific, and (2) calibrating the model parameters. This work provides a framework to address both challenges using a combination of structural imaging at high magnetic fields, stimulus-evoked behavior using a monopolar review, and functional imaging.

### Conductivity Anisotropy

A subject-specific model was created for Monkey L, but the conductivity, fractional anisotropy, diffusion tensors, and conductivity distributions were analyzed for both Monkeys L and P. Fractional anisotropy measured the difference between the three eigenvalues of the diffusion tensor (Pierpaoli et al., 1996). If the diffusion was isotropic (all three eigenvalues are equal), this value became 0. If a large number was calculated, there was high diffusion anisotropy. These values were scaled between 0 (isotropic) and 1 (anisotropic) and displayed as black and white, respectively (**Figure 2C**). In the brainstem region of Monkey L and Monkey P, the fractional anisotropy was found to be highly variable, with values ranging from less than 0.1–0.7. Since the voxel size of the DTI was 1 mm isotropic, each voxel could be composed of multiple fiber tracts, explaining this variability. The highest fractional anisotropy values (∼0.5–0.7 in Monkey L, ∼0.3–0.7 in Monkey P), appeared to correspond to areas of the superior cerebellar peduncle (caudal of the decussation) as well as the medial lemniscus, two of the largest pathways in the brainstem. However, the fractional anisotropy values for the selected slices in **Figure 2C** at the decussation of SCP were small (Monkey L: 0.295, 0.224, Monkey P: 0.181, 0.278, 0.268) and did not vary much from the mean of the surrounding voxels (Monkey L: 0.2867 ± 0.0752, Monkey P: 0.2575 ± 0.0504).

All six parameters of the diffusion tensor are visualized as spherical functions in **Figure 2C**. In both Monkey L and Monkey P, the greatest difference in the overall tensor direction in the brainstem was seen in the area of the ML, where the tensors were primarily dorsal-caudal and were oriented at a 45 degree angle from the surrounding voxels. Additionally, the principal direction (V1) at the decussation of the SCP was oriented medial-lateral, with a 90◦ difference compared to the neighboring voxels. When comparing neighboring voxels elsewhere in the brainstem, many midline voxels displayed at least a 45◦ difference in the longest axis. While some variability between animals was expected, the overall anisotropy in the brainstem was comparable.

Given that the brainstem was composed of a heterogeneous set of nuclei and fiber tracts, we hypothesized that an FEM of the brainstem with anisotropic conductivity would exhibit strong asymmetries in comparison to an otherwise equivalent isotropic model. The conductivity values, σxx, derived from the subjectspecific imaging are shown in **Figure 2A** for several coronal sections throughout the brain. Histograms of the conductivity values calculated from the primary, secondary and tertiary eigenvalues were given for the entire brain, including ventricles, and for the brainstem region around the DBS lead (**Figure 2B**). This region included the pons and part of the midbrain, demarcated as posterior of the substantia nigra. The average calculated conductivity along the main axes in the brainstem was between 0.3235 and 0.4018, just above 0.3 S/m, the value used for the isotropic models.

When conductivity anisotropy was incorporated into the models, the spread of current in the tissue was altered. As seen in **Figure 1D**, the isosurfaces of the electric potential in the isotropic model were spherical, while the isosurfaces in the anisotropic model were non-spherical. This is consistent with previous modeling studies that incorporated anisotropy (Miocinovic et al., 2009). Model predictions for activation threshold were much lower for anisotropic models than for the isotropic model. For example, the threshold for activating 5% of CTG axons using contact 7 was 0.5 mA for the anisotropic model and 1.1 mA for the isotropic model.

The conductivity scaling factor (s) has been previously reported as the range of s = 0.844 ± 0.0545 S·s/mm<sup>3</sup> (Tuch et al., 2001) with varying scaling factors used in other modeling studies (McIntyre et al., 2004b; Butson et al., 2006). To investigate model sensitivity to the conductivity scaling factor within the reported range, activation threshold curves were generated using three values for s (0.79, 0.844, and 0.89). These results are shown in the second column for the ON (**Figure 3**), SCP (**Figure 4**), and CTG (**Figure 5**) tracts assuming an 8 µm diameter axonal fiber for ON and a 2µm diameter for CTG and SCP. Varying the scaling factor ±6.5% resulted in only a minor shift (0.0981 mA) in the threshold for 5% activation of ON fibers using contact 7.

#### Model Parameter Sweep

In addition to investigating model sensitivity to tissue conductivity, other model parameters known to impact the calculation of activation thresholds such as axon diameter (Rattay, 1999) and precise lead location were investigated as well. Previous models of axons in the brainstem region (SCP, medial lemniscus, lateral lemniscus) were modeled with a diameter of 2µm (Zitella et al., 2013). While this is a conservative estimation for SCP and CTG axonal diameters, we also examined the effects on activation thresholds when using 5.7 and 8.7µm axon diameter, which may be more representative of actual axon diameters within the SCP (Hazrati and Parent, 1992) and CTG tracts. The ON tract was also modeled with 2, 5.7, and 8.7µm axon diameter with the latter thought to be the most realistic axon diameter. Atomic force microscopy has shown human oculomotor nerve fibers with much larger diameter fibers, between 10 and 15µm (Melling et al., 2003), which presumably would be slightly smaller in the rhesus macaque. Consistent with the principles of cable models of myelinated axons, the axonal diameter had a large effect on the resultant activation thresholds for all three fiber tracts (**Figures 3**–**5**).

Model sensitivity to the precise position of the DBS lead within the brainstem was investigated by shifting the lead in four directions. Using the same implantation angle as defined by the SWI/CT co-registration process and histological reconstructions from Monkey L, the DBS lead was shifted 0.5 mm anterior, poster, medial, and lateral of the original lead placement (chamber reference). Moving the lead medially increased ON activation, while moving the lead laterally decreased ON activation. At the threshold for contact 7 (1 mA), the model predicted a 27.5%

increase in ON activation for a medial lead location. Overall, anterior and posterior deviation in lead location did not alter the model results for amplitudes near the threshold. Predicted activation from stimulation through contact 5 at 1.4 mA only decreased by 0.7% when the lead was moved in the anterior direction (**Figure 3**).

For contacts 5, 6, and 7, SCP activation increased when the lead position was shifted medial and posterior, while an anterior and lateral lead position decreased activation. SCP activation at 0.9 mA through contact 7 increased the activation by 2.1% (from 0 to 2.1%) when moving the lead posterior. Contact 4 was embedded within SCP, so lead location had minimal effect on SCP activation at lower amplitudes (below 2 mA). However, for amplitudes above 2 mA, the anterior and lateral lead placements decreased activation and the medial lead placement increased activation. For 0.9 mA stimulation through contact 4, change in SCP activation was negligible (∼0.1%) (**Figure 4**).

Due to the anatomy of the CTG, the effect of lead location was different for each active contact. In the posterior direction, there was minimal change in activation for contacts 6 and 7. The same was true for contact 4, at stimulation amplitudes below 1.6 mA. Contact 5 stimulation produced lower activation with a posterior lead placement until stimulation amplitude increased beyond 2 mA, which resulted in a large increase in activation that exceeded the original lead placement results. For all contacts, lateral shift decreased activation and medial shift increased activation, but the magnitude of these changes in activation differed for each contact. For the CTG, moving the lead 0.5 mm medially increased tract activation by 8% when stimulating through Contact 7 at 0.9 mA and 2.6% when stimulating through Contact 4 (**Figure 5**).

# Comparison of ON Model Simulations to Stimulus-Induced Eyelid Flutter

A monopolar review was conducted across all 8 contacts by applying a 20 Hz train of 90µs pulses in increasing amplitude at intervals of 0.1 mA until a motor side effect was observed or the amplitude of stimulation reached 3.5 mA. For the proximal four electrodes (contacts 4–7), a right eyelid flutter was observed at amplitudes at or above 1 mA (**Table 1**) with more proximal contacts requiring higher stimulation amplitudes. At 20 Hz stimulation, the therapeutic PPTg-DBS stimulation frequency, stimulation resulted in an eyelid flutter, while stimulation at higher frequencies (e.g., 130 Hz) resulted in the eyelid remaining elevated. No other overt motor signs were observed at any of the stimulation amplitudes tested for other contacts.

The oculomotor nerve is known to project to the levator palpabrae superioris muscle of the eye, which is responsible for elevation of the upper eyelid (Porter et al., 1989). Multi-compartment axon models were developed for the oculomotor nerve to identify model parameter settings that resulted in the most consistent activation values across the empirical motor threshold amplitude values as defined in **Table 1**. We assumed that the neuron models should

TABLE 3 | Comparison of ON model simulations to behavior thresholds.

predict an activation of 5–15% based on previous models of the corticospinal tract of internal capsule (Chaturvedi et al., 2010) at the experimental motor threshold. Using the ON computational model (diameter = 8.7µm, s = 0.844, original lead location), the percentage of activated axons at motor thresholds was calculated (**Table 3**). For each behavioral threshold, the percent error was calculated as the difference between the experimental threshold and the modelpredicted stimulation amplitude necessary to activate 5% of the axons. There was no error in the model predictions for contact 5, 6, and 7. The percent error for contact 4 was 6.67%.

The neuron modeling results from the anisotropic model resulted in much lower activation thresholds than were predicted from the isotropic model. Moreover, the anisotropic models, in comparison to the isotropic models, resulted in activation thresholds that were more consistent with the thresholds for inducing eyelid flutter (**Figure 3**). For the isotropic models, there was 0% activation of ON at the threshold amplitude for all contacts. Using contact 5, a 5% activation of ON fibers was achieved at 1.25 mA for the anisotropic model, while the isotropic model required 3.4 mA to reach 5% activation. Similarly for contact 6, 1.1 mA activated 8.7% of the axons in the anisotropic model and 2.7 mA was required to activate 8.7% of the axons in the isotropic model. There were equivalent results for the SCP, where the anisotropic model predicted 26% of SCP axons activated at 1.2 mA through contact 4, while the isotropic model predicted 26% of SCP axons activated at 2.5 mA.

# Comparison of Model Simulations to PET Imaging

Two FDG-PET scans, in the context of DBS, were conducted to examine the effects of DBS in the PPTg area, after 0.9 mA stimulation through Contact 7 (configuration 1) and 1.2 mA stimulation through Contact 4 (configuration 2) (**Tables 4, 5**). These were compared to a baseline scan with no stimulation. A sampling of the resultant FDG standard uptake values (FDG-SUV) are shown in **Figure 6**. For the first stimulation configuration, the ventral posteromedial nucleus of thalamus (VPM), which is innervated by CTG (Blumenfeld, 2002), showed an increased FDG-SUV (p = 0.023). Further, descending projections of CTG project to the inferior olivary nuclei (Blumenfeld, 2002), which also showed an increased FDG-SUV

#### TABLE 4 | PET Configuration 1.

in configuration 1 (p = 0.033). This corresponded to an activation of 12.2% of CTG fibers (**Table 3**). For configuration 2, the models predicted no activation of CTG; the PET measured a significant increase in FDG-SUV in the VPM (p = 0.041), but no significant increase in FDG-SUV in the inferior olivary complex (p > 0.05).

Similarly, regions that are innervated by projections from the PPTg showed an increased FDG-SUV, including the centromedian nucleus of thalamus (p = 0.037) during configuration 1. Configuration 2 showed increased FDG-SUV in regions innervated by fiber pathways near PPTg, including the rostral interstitial nucleus of the MLF and the interstitial nucleus of Cajal, which are innervated by MLF. Increased FDG-SUV in downstream targets of the PPTg was also observed in configuration 2, including the globus pallidus, basal amygdala, peripeduncular nucleus, centromedian nucleus, and the STN. The PET results also showed an increase in FDG-SUV in the red nucleus for configuration 2. This was supported by the model predicted activation of the SCP tract for configuration 2, 25.10% (**Table 6**).

# Discussion

Subject-specific computational models are an important tool to better understand the mechanisms of DBS in the brainstem and guide future DBS therapies (Butson et al., 2007; Miocinovic et al., 2009; Chaturvedi et al., 2010; Keane et al., 2012; Lujan et al., 2012; Zitella et al., 2013). In order for these models to be clinically relevant they must provide accurate predictions. While other methods of validation have been applied to

#### TABLE 5 | PET Configuration 2.

FIGURE 6 | PET imaging during PPTg-DBS. (Top) Co-registration of SWI and baseline PET/CT images. The SWI is shown in blue-green cold scale for differentiation from the gray CT. (Left) FDG-SUV during PET configuration 1 (0.9 mA stimulation through

#### TABLE 6 | Model comparison to behavior.

computational models of DBS, no models of the brainstem have yet been rendered subject specific. In this study, we evaluated the sensitivity of a subject-specific model of PPTg-DBS in a nonhuman primate to different model parameters (tissue conductance anisotropy, axonal diameter, and DBS lead location) contact 7), normalized to OFF-DBS. The PET results are overlaid on SWI of Monkey L. (Right) FDG-SUV during PET configuration 2 (1.2 mA stimulation through contact 4), normalized to OFF-DBS.

and compared the results to behavioral and functional imaging measures to determine the most accurate tissue conductance model.

Our previous computational models assumed the DBS lead was surrounded by homogeneous, isotropic tissue with a conductivity of 0.3 S/m (Zitella et al., 2013). Based on the fractional anisotropy results from Monkey L and Monkey P in the brainstem, the mean of the image-based conductivity distribution did deviate from this isotropic conductivity assumption, but was well within an order of magnitude. Since the conductivity scaling factor did not greatly affect the model predictions, the spatial variability of the conductivity (i.e., the distribution of conductivities within the brainstem) proved to have a large effect on the potential distribution around the DBS lead. This

high anisotropy near the lead resulted in lower stimulation amplitudes required to activate nearby axons despite the slightly higher average conductivity in the brainstem. Based on these results, it seems that anisotropy, in conjunction with the average conductivity, plays a role in the ability to activate axons.

Several other computational models of DBS have incorporated anisotropy of tissue conductivity, including models of STN DBS, which assigned typical conductivity values based on the literature (Sotiropoulos and Steinmetz, 2007; Åström et al., 2009). Other studies have converted fractional anisotropy measures from subject DTI to conductivity, including models of DBS in the STN region and in the thalamus (McIntyre et al., 2004a; Miocinovic et al., 2009). Similar to our models, regions with high anisotropy showed greater variability in the voltage isosurfaces and in the activated volume of tissue. However, these studies showed, contrary to our results, that the addition of anisotropy to the model decreased the percentage of axons activated. This difference may relate to axonal fiber orientations relative to the stimulated electrode(s) as well to assumptions of the neuron model parameters.

The PPTg area, similar to other typical DBS target regions, is highly anisotropic. Indeed, the PPTg is surrounded by the spinothalamic tract, CTG, medial lemniscus, lateral lemniscus, and the MLF. The fibers of the SCP are intertwined with the cells of the PPTg, which introduces challenges when attempting to stimulate one pathway over another. The present study showed that the inclusion of anisotropic conductivity is highly important for computational model predictions. This finding suggests that efforts to increase the resolution of fractional anisotropy imaging within the brainstem—through high-field, high angular resolution diffusion imaging (Lenglet et al., 2012), customized head coils (Adriany et al., 2010), and advanced computational reconstruction algorithms (Duarte-Carvajalino et al., 2014) as used in this study—could have significant merit (Novak et al., 2001; Stieltjes et al., 2001; Soria et al., 2011; Aggarwal et al., 2013; Deistung et al., 2013; Gizewski et al., 2014).

The clinical value of computational models lies in their legitimacy, making validation extremely important. Previous studies have confirmed the validity of their model parameters based upon correlations between EMG thresholds and activation of the corticospinal tract of internal capsule (during STN-DBS) (Butson et al., 2006; Chaturvedi et al., 2010), between paresthesia thresholds and activation of the somatosensory representation of thalamus (during Vim-DBS) (Kuncel et al., 2008), and between conjugate eye deviation and oculomotor nerve activation (during STN-DBS) (Butson et al., 2006). Here we extend this approach to the brainstem in the context of DBS targeting the PPTg area. The results showed much better predictions of the activation of the oculomotor nerve axons at stimulation amplitudes necessary to induce eyelid flutter for the anisotropic models. Without being able to measure the magnitude of the eyelid twitch or obtain verbal feedback on the strength of the side effect, this is a positive result. Additional modifications to the model equations, more precise anatomical geometry, and higher resolution DTI could provide more accurate results. The assumption that the conductivity is linearly related to the diffusion tensor eigenvalues may not hold for high resolution (1 mm voxels) DTI within the brainstem and could also be a source of error in the models.

In this study, PET imaging was used as a gross measure of the activation in the area during stimulation to compare the effects of stimulation through different contacts to baseline. PET is a valuable tool that has been used to examine the effects of DBS (Haslinger et al., 2003; Mayberg et al., 2005). The use of PET in the context of PPTg-DBS provided a novel approach to further evaluate the predictive capabilities of the computational neuron models. While the results were consistent, there are several limitations that should be noted. In addition to having only one subject, there was only one scan taken of each configuration (OFF DBS, C4 stimulation, C7 stimulation). Additional small spatial errors could also have been introduced when aligning the INIA19 atlas to the PET/CT. Furthermore, the PET analysis reported here did not account for the precise time from injection to time of scan. This will be incorporated in future studies for more accurate results.

Future studies will also need a larger sample size and expand the model validation methods. Through studies in nonhuman primates, the addition of electrophysiology would provide more insight into the effects of stimulation. The electrophysiological activation thresholds could be compared to the model predictions by recording single-unit spike activity at multiple sites within upstream and downstream targets of fiber pathways coursing near the PPTg.

As DBS techniques continue to advance, new targets are being explored and new lead designs are being developed. There is a growing need for validated computational models to better understand the therapeutic results and titrate stimulation parameters in human patients implanted with DBS systems. This study is the first case of incorporating anisotropic conductivity into subject-specific computational models of DBS in the brainstem. Moreover, the study emphasizes how coupling behavioral metrics and functional imaging data in computational modeling studies can be critical for enhancing the predictive power of the models.

# Acknowledgments

We thank the Center for Magnetic Resonance Research (CMRR) for providing the imaging resources, specifically Bud Grossman of the University of Minnesota Medical School for the use of PMOD. We also thank the Minnesota Supercomputing Institute for providing the computational resources to complete this study. This work was supported by grants from the National Institutes of Health (R01-NS081118, R01-NS085188, P41- EB015894, P30-NS076408, Human Connectome Project U54- MH091657), the Michael J Fox Foundation, and the National Science Foundation (IGERT DGE-1069104 to LZ and GRFP 00006595 to BT).

# References

treatment of intractable proximal tremor. J. Neurosurg. 99, 708–715. doi: 10.3171/jns.2<IP_ADDRESS>8

**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.

Copyright © 2015 Zitella, Teplitzky, Yager, Hudson, Brintz, Duchin, Harel, Vitek, Baker and Johnson. 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.

# Motor planning under temporal uncertainty is suboptimal when the gain function is asymmetric

Keiji Ota1, 2 \*, Masahiro Shinya<sup>1</sup> and Kazutoshi Kudo<sup>1</sup> \*

*<sup>1</sup> Laboratory of Sports Sciences, Department of Life Sciences, Graduate School of Arts and Sciences, The University of Tokyo, Tokyo, Japan, <sup>2</sup> Research Fellow of Japan Society for the Promotion of Science, Tokyo, Japan*

For optimal action planning, the gain/loss associated with actions and the variability in motor output should both be considered. A number of studies make conflicting claims about the optimality of human action planning but cannot be reconciled due to their use of different movements and gain/loss functions. The disagreement is possibly because of differences in the experimental design and differences in the energetic cost of participant motor effort. We used a coincident timing task, which requires decision making with constant energetic cost, to test the optimality of participant's timing strategies under four configurations of the gain function. We compared participant strategies to an optimal timing strategy calculated from a Bayesian model that maximizes the expected gain. We found suboptimal timing strategies under two configurations of the gain function characterized by asymmetry, in which higher gain is associated with higher risk of zero gain. Participants showed a risk-seeking strategy by responding closer than optimal to the time of onset/offset of zero gain. Meanwhile, there was good agreement of the model with actual performance under two configurations of the gain function characterized by symmetry. Our findings show that human ability to make decisions that must reflect uncertainty in one's own motor output has limits that depend on the configuration of the gain function.

#### Edited by:

*Vincent C. K. Cheung, The Chinese University of Hong Kong, Hong Kong*

#### Reviewed by:

*Athena Akrami, Princeton University - Howard Hughes Medical Institute, USA Shih-Wei Wu, National Yang-Ming University, Taiwan*

#### \*Correspondence:

*Keiji Ota and Kazutoshi Kudo, Laboratory of Sports Sciences, Department of Life Sciences, Graduate School of Arts and Sciences, The University of Tokyo, Komaba 3-8-1, Meguro, Tokyo 153-8902, Japan<EMAIL_ADDRESS><EMAIL_ADDRESS>
> Received: *09 February 2015* Accepted: *25 June 2015* Published: *15 July 2015*

#### Citation:

*Ota K, Shinya M and Kudo K (2015) Motor planning under temporal uncertainty is suboptimal when the gain function is asymmetric. Front. Comput. Neurosci. 9:88. doi: 10.3389/fncom.2015.00088* Keywords: decision making, risk-sensitivity, Bayesian decision model, response variance, coincident timing task

# Introduction

In highly skilled movement, especially in sports, decision making is important for superior performance. For example, a tennis player requires a spatial action plan about where in a court they should aim; a ski jumper requires a temporal action plan about when they should take off. An executed action is associated with a gain/loss. In ski jumping, the take-off jump should be as close to the edge of the ramp as possible to get the best jump length, while take-off too early or too late decreases jump length (Müller, 2009). However, in whatever action they plan, an executed action is not always equal to the planned one because of motor variability (Schmidt et al., 1979; Kudo et al., 2000; van Beers et al., 2004). Thus, both gain/loss associated with action and uncertainty in motor output should be considered for better decision making.

The mathematical method for selecting an optimal plan under conditions of limited uncertainty is known as statistical decision theory (Berger, 1985; Maloney and Zhang, 2010). In particular, Bayesian decision theory, which is a part of statistical decision theory, is a method for optimizing the expected gain/loss.The expected gain/loss is calculated by integrating the gain/loss function assigned to a certain action over a probability distribution of an executed action given a planned action. The Bayesian decision maker plans the action that optimizes the expected gain/loss for any combination of gain/loss function and motor variability (Hudson et al., 2012).

Previous motor control studies have evaluated the optimality of human motor decision making by comparing Bayesian ideal performance with actual human performance (Trommershäuser et al., 2003a,b, 2005; Wu et al., 2006; Hudson et al., 2012; O'Brien and Ahmed, 2013). Some reports have concluded that humans can plan actions that are optimal when considering their own motor variability (Trommershäuser et al., 2003a,b, 2005; Hudson et al., 2012), while other reports have concluded that humans cannot compute the movement strategy that maximizes the expected gain in the presence of such variability (Wu et al., 2006; O'Brien and Ahmed, 2013). Thus, there is inconsistency in published claims concerning the optimality of human action planning.

Two possible factors,—differences in experimental design, and differences in energetic cost—could account for this inconsistency. First, previous studies have reported optimal or suboptimal action plans with different movements and different configurations of the gain function. For example, Trommershäuser et al. (2003a,b, 2005) and Hudson et al. (2012) have demonstrated optimality in a pointing plan under a gain function in which the magnitude of gain/loss remains constant. In contrast, O'Brien and Ahmed (2013) have shown suboptimal reaching and whole-body movement strategies under an asymmetric gain function in which seeking higher values of gain brings participants closer to scoring zero gain ("falling over the cliff "). Therefore, we cannot directly evaluate the relationship between the optimality of the action plan and the configuration of the gain function because the experimental designs among previous studies differed.

Second, previous studies have mainly treated pointing or reaching movement as executed action (Trommershäuser et al., 2003a,b, 2005; Wu et al., 2006; Hudson et al., 2012; O'Brien and Ahmed, 2013). In reaching and pointing movements, energetic cost is proportionally larger as the distance of the required movement is made longer. Because large energetic cost requires participant large effort, energetic cost could be a factor disturbing the measured optimality of action strategies. For example, Hudson et al. (2012) have reported that a discrepancy between ideal and actual performance emerged when optimal but largecost movements were required during obstacle avoidance (i.e., large excursions).

Here, we used a coincident timing task requiring decision making and compared the Bayesian ideal performance with actual human performance under four different configurations of gain function including those used O'Brien and Ahmed (2013). In the coincident timing task, energetic cost is constant because the participant just presses a button whatever strategy he/she selects. Thus, we can directly evaluate the relationship between the optimality of action plans and the configuration of the gain function excluding the factor of energetic cost as a possible reason for any discrepancy found. We, in fact, found good agreement between the ideal timing strategy and the actual strategy under a symmetric configuration. However, a discrepancy was found under asymmetric configurations. We will discuss possible explanations for this discrepancy. We also observed that larger trial-by-trial compensation occurred following miss trials than after success trials even though the experienced response errors were of the same magnitude.

# Methods

# Participants

Thirty-seven healthy right-handed adults participated in the experiment. Sixteen participants (10 male, 6 female; mean age 28.1 ± 7.8 years) performed Experiment 1, twelve (10 male, 2 female; mean age 22.8 ± 2.8 years) participants performed Experiment 2, and the remaining nine participants (7 male, 2 female; mean age 21.3 ± 2.2 years) performed Experiment 3. All participants were unaware of the purpose of the experiment. This study was approved by the Ethics Committee of the Graduate School of Arts and Sciences, the University of Tokyo.

# Experimental Task

**Figure 1** shows the time sequence of our basic experimental task. First, a warning tone was presented to ready the participants for an upcoming trial. Then, a visual cue was presented on a computer screen as a starting signal (14 inches, 1600×900 pixels, refresh frequency 60 Hz; Latitude E5420, Dell, Round Rock, TX, USA). The participants were instructed to press a button after presentation of the visual cue. The response time was recorded as the button-press time relative to the onset time of the visual cue. In each trial, the participants gained a point based on "gain function," a function that translated the response time to a certain number of points. The details of the gain function are explained in the following section. The foreperiod (interval between the warning tone and the visual cue) was randomly varied between 800 and 1200 ms in steps of 100 ms. The target time was 2300 ms after the visual cue and was fixed throughout the experiment. In our experiment, the target time was associated with gaining 100 points but was not necessarily the time when the participants should respond (see below). The inter-trial interval was 2000 ms. All computerized events were controlled by a program written with LabVIEW software (National Instruments, 2011 Service Pack 1, Austin, TX, USA).

## Experimental Condition and Procedure

In each experimental condition, the participants were required to make a decision about when to press a button to maximize the total gain in 100 trials. The gain for a trial was a function of response time, termed the "gain function." Four conditions were tested, corresponding to different gain functions. The first was characterized as the No Risk condition, which employed a symmetric gain function (**Figure 2A**). In this condition, a gain for a trial (G) was determined from the following equation.

$$G(t) = \begin{cases} \begin{array}{c} \frac{1}{23}t, \text{ if } t \le 2300 \\\ -\frac{1}{23}t + 200, \text{ if } t > 2300 \end{array} \end{cases} \tag{1}$$

In the above equation, t represents a response time in milliseconds. When the participants responded earlier than the

target time, they received a number of points that was a positive linear function of response interval. When the participants responded later than the target time, they received a number of points that was a negative linear function of response interval. A maximum-possible one-trial gain of 100 points could be obtained by responding exactly at the target time.

The second condition was the Step condition, which also had a symmetric gain function (**Figure 2B**). In the Step condition, a period of constant gain was on the edge of risk of zero gain both at its start and at its termination as represented by the following equation.

$$G(t) = \begin{cases} \begin{array}{c} 0, \text{ } \acute{y} \\ 100, \text{ } \acute{y} \\ 0, \text{ } \acute{y} \end{array} 1900 \stackrel{\le}{\le} t \stackrel{\le}{\le} 2700 \\\ \text{ } \text{ } \acute{y} \end{array} \tag{2}$$

The participants received 100 points if they responded within ±400 ms of the target time. However, zero points were given if they responded at less than target −400 ms or at more than target +400 ms. We termed these eventualities "miss trials" and when they occurred, the participants were cautioned by an unpleasant alarm and a flashing red lamp on the screen. The volume of the alarm was 71.5 ± 0.4 dB.

The third condition was characterized as the Riskafter condition, which employed an asymmetric gain function (**Figure 2C**). In the Riskafter condition, single-trial gain rose linearly as the target time approached, then plunged to zero after the target time and then remained zero thereafter as represented by the following equation.

$$G(t) = \begin{cases} \frac{1}{23}t, & \text{if } t \le 2300\\ 0, & \text{if } t > 2300 \end{cases} \tag{3}$$

Earlier than the target time, the gain function looks the same as that in the No Risk condition. However, zero points were given if the participants responded after the target time. If they missed, they received the same penalty as in the Step condition. This configuration has been used in O'Brien and Ahmed (2013).

The last condition was the Riskbefore condition in which the gain function of the Riskafter condition is mirror-imaged across the target time (**Figure 2D**). The gain for a trial was determined from the following equation.

$$G(t) = \begin{cases} 0, \text{ if } t < 2300 \\ -\frac{1}{23}t + 200, \text{ if } t \ge 2300 \end{cases} \tag{4}$$

In contrast to the Riskafter condition, zero points were given if the participants responded before the target time in the Riskbefore condition. Again, if they missed, they received the same penalty as in the Step condition. For trial-by-trial compensation analysis described below, we defined miss trials as trials when the participants received zero points, success trials as any other trials in all four conditions. Of note, all the trials resulted in success in the No Risk condition because the range in success trials was enough large (i.e., it was 4600 ms).

In each trial, we provided the participants with feedback information consisting of the relative response time calculated by response time–the target time, the gain for the trial, and the cumulative total gain. All the participants performed 10 trials for practice. This practice session was conducted to give them a feel for the length of time from the visual cue to the target time. In this session, we only provided them with the relative response time (i.e., no gain function was applied). After the practice session, all the participants performed 100 trials in the No Risk condition as a first experimental session. In the second experimental session, the participants who were assigned to Experiment 1 performed 100 trials in the Riskafter condition. The participants who were assigned to Experiment 2 performed 100 trials in the Riskbefore condition. Those who were assigned to Experiment 3 performed 100 trials in the Step condition. Each condition (No Risk and Riskbefore/Riskbefore/Step) was conducted as separate experimental session. The participants rested for several tens of seconds between sessions.

Before running each session, we explained the structure of the gain function with a figure visualizing it (see Supplementary Figure 1). Thus, the structure of the gain function was known to the participants before testing. In the figure, following information was also included: 100 points could be gained when the relative response time was 0 ms in the No Risk, Riskafter, and the Riskbefore conditions, and when the relative response time was within ±400 ms of the target in the Step condition. However, the length of time from the visual cue to the target time was not described; thus, the participants did not know that it was 2300 ms. Also, before performing the No Risk condition, the participants did not know the gain structure that would be used in the next condition.

We instructed them to maximize total gain in each condition. Thus, they were required to make a decision about when to press a button to maximize the total gain. Actual monetary rewards/penalties were not used (see limitation related to this experimental procedure).

#### Model Assumptions

We calculated the ideal strategy that maximizes the expected gain by a Bayesian decision-theoretic approach for each participant and for all conditions (Hudson et al., 2012). Our model consisted of two sets and two functions. The two sets were: possible response strategy T (motor decision), and executed response time t (motor output). The two functions were: probability distribution of executed response P (t|T), and gain function G(t). Given a particular planned strategy, a particular response is stochastically executed. This is considered the uncertainty in motor output. In this study, we assumed that the produced response time t is distributed around the planned response time T according a Gaussian distribution (see Supplementary Table 1) as follows.

$$P\left(t|T\right) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{\left(t-T\right)^2}{2\sigma^2}\right] \tag{5}$$

Then, given execution of a particular response time, the gain is given according to the gain function G(t). Given both P (t|T) and G(t), the expected gain EG(T) as a function of planned response time T is calculated by the following equation.

$$EG(T) = \int\_{-\infty}^{\infty} G(t) \cdot P(t|T)dt\tag{6}$$

Once we had measured the response variance σ for each participant and condition, we could calculate the optimal mean response time T <sup>∗</sup> by maximizing Equation (6). A Bayesian decision maker chooses a response time T ∗ for any given gain function G (t) and response variance σ. This can be regarded as a theoretical risk-neutral optimal response.

# Estimation of 95% Confidence Interval of Optimal Response Time

Furthermore, we estimated the 95% confidence interval of the optimal mean response time T <sup>∗</sup> by bootstrapping (3000 resamples) for the Riskafter and the Riskbefore conditions. We then examined whether the actual response time is within this 95% confidence interval in the Riskafter condition and in the Riskbefore condition. In a range from σ = 0 to 0.4 in steps of σ = 0.001, we first calculated each optimal mean response time Xopt1(σ0), Xopt2(σ0.001), · · ·, Xopt400(σ0.4) by maximizing Equation (6). Here focusing on Xopt1(σ0), we simulated 100 trials of a task execution responding by this optimal mean response time and having this response variance (in this case Xopt1 = 2300, σ = 0) using a MATLAB randn function.  For each bootstrap sample x<sup>b</sup> (b = 1, 2, · · ·, 3000), we calculated the average value of its samples <sup>b</sup>µ<sup>b</sup> <sup>=</sup> 1 <sup>100</sup> · P<sup>100</sup> i=1 (xbti). After sorting these average samples bµb (b = 1, 2, · · ·, 3000) in ascending order, we defined a 2.5% and a 97.5% point in these samples as the 95% confidence interval in optimal mean response time T ∗ . By repeating the above processes from Xopt<sup>1</sup> (σ0) to Xopt400(σ0.4), we estimated the 95% confidence intervals of each optimal mean response time.

If the observed mean response times were within these 95% confidence intervals, we would conclude that the participant plans optimal and risk-neutral timing strategies. In the Riskafter condition, observed times longer than the confidence intervals would indicate suboptimal and risk-seeking strategies. Observed times shorter than the confidence intervals would indicate suboptimal and risk-averse strategies.

# Optimal Response Times Calculated from the Measured Distributions

Although we had confirmed that the response distributions were Gaussian (see Supplementary Table 1), we also calculated optimal mean response time for the Riskafter and the Riskbefore conditions using the measured response distributions. In the Riskafter condition (**Figure 3A**, upper panel), once we had obtained the response distribution we simply shifted it back and forth to identify the maximal total gain for that distribution (**Figure 3A**, lower panel). In **Figure 3A**, we show the case of shifting the measured distribution back. We defined the optimal mean response time as the mean response time of the optimized distribution (gray solid line in **Figure 3A**). The estimated optimal mean response time was always earlier than the target time in the Riskafter condition. The difference between the estimated optimal response time and the target time reflected each participant's own variance in response time, (i.e., the larger one's variance, the

earlier the optimal response time, and vice-versa. This effect is visualized in **Figure 5** as the solid curve).

In **Figure 3B**, the black arrow indicates the total gain when the measured distribution is not shifted and the gray arrow indicates the highest total gain possible for that distribution under shifting. The two-headed arrows in **Figures 3A,B** represent the optimal shift size.

Finally, we compared the estimated optimal mean response times with the observed mean response times (black solid line in **Figure 3A**). In this case, the observed time was closer to the target time than optimal, indicating risk-seeking. We also applied this approach to the Riskbefore condition. We compared the mean, as opposed to median optimal and observed times because the measured distributions were Gaussian.

## Trial-By-Trial Compensation Strategy

In addition to determining the response time strategies based on all trials, we examined compensation against the most recent response error based on a trial-by-trial analysis. These results were then compared between/within the No Risk and the Riskafter/Riskbefore conditions. The magnitude of the response error experienced on the current trial is known to influence the response in the following trial (Thoroughman and Shadmehr, 2000; Scheidt et al., 2001). Thus, the compensation size in the following trial can be proportional to the current magnitude of response error. Additionally, it has been shown that humans adjust future motor behavior according to rewarding and nonrewarding outcomes experienced (Wrase et al., 2007). Therefore, in addition to the compensation strategy against response error, we hypothesized that the compensation size following miss trials would be larger in the Riskafter/Riskbefore conditions than the compensation size following success trials in the No Risk condition.

We defined the compensation on the current trial, trial n, by subtracting the response time on the current trial, from that on the following trial, trial n+1, as in Equation (7).

$$Compensation\_n = RT\_{n+1} - RT\_n \tag{7}$$

We supposed that the compensation occurs around mean response time in both conditions, thus we defined response error as response time–mean response time in this analysis. The compensation size was anticipated to depend on the magnitude of response error (in other words, the magnitude of deviation between the current response time and mean response time).