# THALAMIC FUNCTION - BEYOND A SIMPLE RELAY

EDITED BY: Vincenzo Crunelli, William Martin Connelly and W. Martin Usrey PUBLISHED IN: Frontiers in Neural Circuits

#### *Frontiers Copyright Statement*

*© Copyright 2007-2016 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

> *All copyright, and all rights therein, are protected by national and international copyright laws.*

*The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88919-842-9 DOI 10.3389/978-2-88919-842-9

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# **THALAMIC FUNCTION - BEYOND A SIMPLE RELAY**

Topic Editors: **Vincenzo Crunelli,** Cardiff University, UK

**William Martin Connelly,** Australian National University, Australia **W. Martin Usrey,** University of California, Davis, USA

The thalamus is often described as a relay. Typified by sensory pathways, this concept leads to thalamic nuclei being viewed as areas that passively streams information from a single source to the cortex, without affecting the nature of that information. However, diverse intrathalamic connections, the varying synaptic and membrane properties of thalamic neurons and the large number of inputs from non-sensory sources make the idea that the thalamus is just a passive relay unlikely. Furthermore, a large number of thalamic nuclei are not primarily driven by sensory signals nor do they exclusively target the cortex, meaning the thalamus must do more than simply pass sensory signals to the cortex. Finally, there is a wealth of research demonstrating that the thalamus does indeed function in ways that are not captured by the concept of a simple relay. So why, given all of this, is the primary paradigm for describing the thalamus, a relay? This Research Topic covers original research, reviews and hypotheses on thalamic function that explore the concept that the thalamus performs computational tasks other than simply passively relaying information.

**Citation:** Crunelli, V., Connelly, W. M., Usrey, W. M., eds. (2016). Thalamic Function - Beyond a Simple Relay. Lausanne: Frontiers Media. doi: 10.3389/978-2-88919-842-9

# Table of Contents

*153 The Slow Oscillation in Cortical and Thalamic Networks: Mechanisms and Functions*

Garrett T. Neske

François David, Vincenzo Crunelli, Nathalie Leresche and Régis C. Lambert

*218 On Parallel Streams through the Mouse Dorsal Lateral Geniculate Nucleus* Daniel J. Denman and Diego Contreras

# The olfactory thalamus: unanswered questions about the role of the mediodorsal thalamic nucleus in olfaction

### Emmanuelle Courtiol 1,2\* and Donald A. Wilson1,2

1 Emotional Brain Institute, Nathan Kline Institute for Psychiatric Research, Orangeburg, NY, USA, 2 Department of Child and Adolescent Psychiatry, New York University Langone Medical Center, NY, USA

The mediodorsal thalamic nucleus (MDT) is a higher order thalamic nucleus and its role in cognition is increasingly well established. Interestingly, components of the MDT also have a somewhat unique sensory function as they link primary olfactory cortex to orbitofrontal associative cortex. In fact, anatomical evidence firmly demonstrates that the MDT receives direct input from primary olfactory areas including the piriform cortex and has dense reciprocal connections with the orbitofrontal cortex. The functions of this olfactory pathway have been poorly explored but lesion, imaging, and electrophysiological studies suggest that these connections may be involved in olfactory processing including odor perception, discrimination, learning, and attention. However, many important questions regarding the MDT and olfaction remain unanswered. Our goal here is not only to briefly review the existing literature but also to highlight some of the remaining questions that need to be answered to better define the role(s) of the MDT in olfactory processing.

#### Edited by:

William Martin Connelly, Cardiff University, UK

#### Reviewed by:

Jane Plailly, Centre National de la Recherche Scientifique, France Burton Slotnick, American University, USA

#### \*Correspondence:

Emmanuelle Courtiol, Emotional Brain Institute, Nathan Kline Institute for Psychiatric Research, 140 Old Orangeburg Road, Orangeburg, NY 10962, USA<EMAIL_ADDRESS>
Received: 16 June 2015 Accepted: 31 August 2015 Published: 18 September 2015

#### Citation:

Courtiol E and Wilson DA (2015) The olfactory thalamus: unanswered questions about the role of the mediodorsal thalamic nucleus in olfaction. Front. Neural Circuits 9:49. doi: 10.3389/fncir.2015.00049 Keywords: olfaction, mediodorsal thalamus, dorsomedial thalamus, piriform cortex, odor response

### Introduction

The thalamus is a crucial crossroad structure in the brain that is recognized as a major contributor to the following functions: sensory perception, attention, sleep and arousal, memory, and cognition. Thalamic nuclei can be divided into (at least) two categories: first-order and higher order thalamic relays (Guillery, 1995). The first category, sensory recipient thalamic relays, processes information arriving from the periphery. The second category, higher order thalamic relays, processes information sent from many cortical areas. Higher order thalamic relays are key structures in cortico-thalamo-cortical networks (Sherman and Guillery, 2002; Mitchell et al., 2014; Saalmann, 2014).

The mediodorsal thalamic nucleus (MDT) is an example of a higher-order thalamic relay (Mitchell and Chakraborty, 2013). The MDT receives inputs from a wide variety of brain areas including cortical structures (notably the prefrontal cortex), brainstem structures, basal forebrain structures, and other thalamic nuclei (Groenewegen, 1988; Kuroda and Price, 1991a,b; Ray and Price, 1992; Guillery, 1995; Kuroda, 1998). In return, the MDT projects massively to the prefrontal cortex (Leonard, 1969; Krettek and Price, 1977a). The cytoarchitecture and the topographical distribution of the different inputs and outputs have led to the separation of the MDT into three subnuclei in the rat—medial, central, and lateral (Krettek and Price, 1977a; Groenewegen, 1988). The dense reciprocal connections between the MDT and the prefrontal cortex have placed the MDT as a critical structure in the study of cognitive processes. In fact, lesions of the MDT in monkeys and rats are associated with a wide range of cognitive deficits: mnesic deficits, deficits in stimulus-outcome associations, deficits in representation of outcome value, and deficits in action-outcome association (Corbit et al., 2003; Mitchell and Gaffan, 2008; Ostlund and Balleine, 2008; Baxter, 2013; Mitchell and Chakraborty, 2013; Alcaraz et al., 2014; Mair et al., 2015). Electrophysiological recordings of the MDT have also demonstrated the contribution of the MDT in working memory, behavioral flexibility, goaldirected behavior, and stimulus reward-association (Oyoshi et al., 1996; Kawagoe et al., 2007; Yu et al., 2012; Han et al., 2013; Parnaudeau et al., 2013; Mair et al., 2015). The role of the MDT in cognition is thus increasingly well established (reviewed in Baxter, 2013; Funahashi, 2013; Mitchell and Chakraborty, 2013; Mitchell et al., 2014; Mitchell, 2015). In addition to these cognitive functions, the MDT also has a sensory component as the olfactory thalamus. As described below, there are firm anatomical and physiological data demonstrating the relationships among the olfactory cortex, the MDT, and the orbitofrontal associative cortex. These connections are particularly intriguing as they bring together one of the most phylogenetically oldest sensory systems with one of the more recently evolved cortical structures.

### Anatomy of the Olfactory Thalamus

The thalamus is the major source of sensory information to the primary sensory cortex for all of the senses except olfaction. In fact, olfactory sensory neurons send their axons directly to the olfactory bulb which in turn projects to the primary olfactory cortex—a region including the piriform cortex, the anterior olfactory nucleus, the olfactory tubercle, the cortical nucleus of the amygdala, and the lateral entorhinal cortex (Price and Powell, 1971; Haberly and Price, 1977; **Figure 1A**). While there is no direct input from the olfactory sensory neurons to the thalamus, the MDT both receives and sends information to primary as well as secondary olfactory areas. An example of a major secondary olfactory area is the orbitofrontal cortex which has strong reciprocal connections with both the MDT and piriform cortex (Illig, 2005). While this review focuses on the MDT, the submedial nucleus of the thalamus also receives olfactory inputs (Price and Slotnick, 1983; Price, 1985).

### Olfactory Afferents

Powell et al. (1963) was one of the first to reveal the relationship between the olfactory pathway and the MDT by showing axonal fiber degeneration in the MDT following lesions in the piriform cortex. In addition to the piriform cortex, the MDT also receives direct input from the olfactory tubercle, the basolateral and cortical nuclei of the amygdala, the lateral entorhinal cortex, the anterior olfactory nucleus, the endopiriform nucleus, and the orbitofrontal cortex (**Figure 1A**). The MDT is thus the target of all the primary olfactory areas (e.g., piriform cortex) as well as some secondary olfactory areas (e.g., orbitofrontal cortex). Of note, the olfactory projections are topographical and involve two distinct subregions of the MDT: the medial and central subnuclei [in rats: (Powell et al., 1963; Heimer, 1968; Krettek and Price, 1974, 1977b; Inagaki et al., 1983; Price and Slotnick, 1983; Price, 1985; Cornwall and Phillipson, 1988; Kuroda and Price, 1991a,b; Kuroda et al., 1992a,b; Kowianski et al., 1999; Bay and Cavdar, 2013; Wilson et al., 2014); in hamsters: (Ferrer, 1969); **Figure 1A**]. The posterior piriform cortex, the anterior olfactory nucleus, the basolateral and cortical nuclei of the amygdala, the agranular insular areas, and the lateral entorhinal cortex project more medially in the MDT while the rostral piriform cortex (deep layers), the ventral and lateral orbital areas, the olfactory tubercle (polymorphic area), and the ventral endopiriform nucleus project mainly to the central region of the MDT [(Krettek and Price, 1974; Inagaki et al., 1983; Price and Slotnick, 1983; Price, 1985; Ray and Price, 1992; Bay and Cavdar, 2013); a small number of neurons from the endopiriform nucleus also project to the medial MDT (Cornwall and Phillipson, 1988)]. The detailed synaptic organization of the olfactory projections to the MDT is still poorly known. Kuroda et al. (1992a,b) described two types of axon terminals (large and small presynaptic terminals) from the piriform cortex to the central MDT which both appear to be excitatory. However, those two types of axon terminals still need to be characterized physiologically to identify whether or not the cells of origin differ. The targeting of specific cell types with cell-specific viral manipulations (e.g., optogenetics) may help to answer this question.

### Olfactory-Related Efferents

The MDT is known to be the origin of dense projections to the frontal cortex in the rat (Groenewegen, 1988). Furthermore, topographical projections from the MDT to olfactory-related structures, including the orbitofrontal cortex and the amygdala, have been established (Krettek and Price, 1977a; Ray and Price, 1992; **Figure 1A**). Interestingly, Kowianski et al. (1999) demonstrated that the endopiriform nucleus receives input from the MDT, though this has not been reported elsewhere.

In essence, the medial subnucleus of the MDT projects to prelimbic and dorsal agranular insular areas, as well as to the basolateral amygdala. The central subnucleus projects to the lateral part of the orbitofrontal cortex and the ventral part of the agranular insular area. Finally, the lateral subnucleus projects to the anterior part of the cingular area, the medial precentral area, and is reciprocally connected with orbital areas (Krettek and Price, 1977a; Groenewegen, 1988; Ray and Price, 1992). The projections of the MDT to the orbitofrontal cortex and basolateral amygdala are of great interest as those two areas are strongly involved in olfactory perception and odor-guided behavior [in rats: (Schoenbaum et al., 1999; Sevelinges et al., 2004; Feierstein et al., 2006; Roesch et al., 2007; Chapuis et al., 2009); in monkeys: (Tanabe et al., 1975; Rolls et al., 1996); in humans: (Jones-Gotman and Zatorre, 1993; Zald and Pardo, 1997; Gottfried and Zelano, 2011)].

The different anatomical studies provide strong evidence establishing the relationship between the olfactory pathway and the MDT. However questions related to the cells of origin and ultrastructural and synaptic organizations of the olfactory afferents in the MDT, as well as the neurotransmitters involved,

still need to be investigated. Furthermore, given the recent demonstration that piriform cortical neurons projecting to orbitofrontal cortex may be non-randomly spatially organized, more detailed analysis of the olfactory cortex-MDT projection is warranted (Chen et al., 2014).

### Electrophysiological Studies of the Olfactory Thalamus

As a first step in understanding the contribution of the MDT in olfactory processing, it is important to characterize how olfactory information is encoded in the MDT. Here, we will describe the physiological responses of the MDT to olfactory stimulation. These data are based on responses recorded in the central and medial portions of the MDT.

First, evoked potentials and extracellular unitary responses in the MDT following the electrical stimulation of the olfactory bulb or lateral olfactory tract of various species have been described [in rats, central subnucleus: (Price and Slotnick, 1983; Price, 1985); in monkeys, medial subnucleus: (Benjamin and Jackson, 1974; Yarita et al., 1980; Takagi, 1986); in rabbits, medial subnucleus: (Jackson and Benjamin, 1974; Imamura et al., 1984)], with percentages of olfactory regiondriven MDT units ranging from 16% to approximately 70% [(Imamura et al., 1984): 87/538 units in rabbits and Benjamin and Jackson (1974): approximately 180 units recorded in monkeys and 127 responsive units were localized in medial MDT].

Second, MDT units can respond to various odorant categories including biological, monomolecular, and mixture odorants [in rats, we observed odor-responsive units in medial, central and lateral MDT, but the boundary delimitation of the different subnuclei was not always clear and may account for the observation of odor-responsive units in the lateral MDT: (Courtiol and Wilson, 2014); in rabbits (Imamura et al., 1984); in monkeys: (Yarita et al., 1980; Takagi, 1986); in cats: (Motokizawa, 1974); **Figure 1B1**]. We observed that 51% of rat MDT units were odor-responsive, Motokizawa (1974) reported 44% of odor-responsive units in cats and Imamura et al. (1984) reported 55% of odor-responsive units in rabbits, although this last percentage was calculated only on MDT neurons responding to lateral olfactory tract shocks (48/87). The tuning of MDT neurons to odorants seems to be dependent on the species, odorant set, and level of anesthesia used. In fact, Courtiol and Wilson (2014) reported that, in anesthetized rats, 63% of units responded to only one odor out of the odor set used. Yet, this tuning has been shown to be more broad in awake monkey and anesthetized rabbits, with 41.5% of MDT neurons responding to four odors and more than 80% responding to two odors or more, respectively (Yarita et al., 1980; Imamura et al., 1984). Nevertheless, two common features of these reports are that MDT units: (1) primarily display excitatory activities in response to odors although suppressive responses can also be observed (**Figure 1B1**) and (2) are similarly responsive to biological, monomolecular or mixture odorants. In addition, our work and the work of others reported that some MDT units can display respiration-linked activity (Imamura et al., 1984; Courtiol and Wilson, 2014). Interestingly, at least in anesthetized rabbits and unanesthetized monkeys, the majority of odor responsive neurons are sensory specific as they do not respond to visual, auditory nor somatic stimulation (Yarita et al., 1980; Imamura et al., 1984).

Finally, the response of the MDT to odorant stimulation can also be recorded at the network level. In recording the local field potentials in the MDT of urethane-anesthetized rats (Courtiol and Wilson, 2014; **Figure 1B2**), we observed that odorants induce the conjoint emergence of beta frequency oscillations in the MDT and piriform cortex. Interestingly, a subset of MDT units fire in phase with the beta frequency oscillations recorded in the piriform cortex. These beta oscillations may offer an effective mechanism for olfactory information transmission between the piriform cortex and the MDT (Tallon-Baudry et al., 2001).

Taken together, these studies reveal that the MDT can respond to and encode odorant information in a manner similar to other primary and secondary olfactory structures. However, future studies will need to determine: (1) the contribution of each of the olfactory inputs to the MDT response; (2) the impact of the MDT on downstream targets such as the orbitofrontal cortex; and (3) given the rich variety of non-olfactory inputs to MDT, how MDT neurons contribute to multi-sensory associations and contextual effects on odor perception. Performing multi-site unitary recordings and using a large set of odorants in behaving animals may achieve this.

### The Still Unclear Role of the MDT in Olfaction

Beyond its basic odor responsiveness, the role of the MDT in olfaction remains unclear. In this last section, we will review the different studies involving the MDT in olfaction, point to some common threads among the available literature, and highlight remaining questions.

Studies of the effect of damage to the MDT in both humans and animal models have provided some useful information about its role in olfaction (Tham et al., 2009). The results of these studies have demonstrated that both humans and animal models with MDT damage are not anosmic and do not present deficits in olfactory detection [in rats and hamsters: (Eichenbaum et al., 1980; Sapolsky and Eichenbaum, 1980); in humans: (Potter and Butters, 1980; Sela et al., 2009)]. Furthermore, humans with Korsakoff 's disease presenting MDT damage are not anosmic but odor detection effects vary between studies, probably depending on the extent of the damage (Jones et al., 1975; Potter and Butters, 1980; Pol et al., 2002). While the results of those studies have shown that MDT damage does not affect olfactory detection, they support the fact that MDT lesions do affect other olfactory functions including olfactory perception, discrimination, learning, and attention. In respect to olfactory perception, Sapolsky and Eichenbaum (1980) demonstrated that MDT-lesioned hamsters show distorted odor preference, i.e., less interest in female and male odors and reduced preference for genital sniffing, leading to maladaptive sexual behaviors. Altered odor preference was also reported in humans. In fact, patients with damage in the MDT present altered olfactory hedonic perception (Rousseaux et al., 1996; Asai et al., 2008; Sela et al., 2009). These results are interesting given the reciprocal connections between the MDT and the amygdala, as the neuronal activity in the amygdala has been shown to be directly influenced by the hedonic valence of olfactory stimuli in humans (Zald and Pardo, 1997). The MDT may thus be part of a network, including the amygdala, involved in the coding of the hedonic valence of the olfactory stimuli.

With respect to discrimination, Eichenbaum et al. (1980) showed that rats with lesions of the MDT exhibit deficits in difficult odor discriminations. For instance, rats with MDT lesions need more trials to reach the discrimination criterion when the task difficulty is increased by using either novel stimuli or perceptually similar stimuli (Eichenbaum et al., 1980; Slotnick and Risser, 1990). However, these deficits can be temporary and alleviated if the animals receive intensive training (Staubli et al., 1987). Deficits in odor discrimination, as well as in odor identification, were also reported in humans with thalamic damage (Sela et al., 2009; Tham et al., 2011a). Notably, Tham et al. (2011b) used both olfactory and visual discrimination tests and demonstrated that patients with MDT damage presented impaired performance selectively for olfactory discrimination compared to visual discrimination. The effects of MDT lesions extend beyond olfactory discrimination given that deficits in olfactory learning have also been observed. For example, rats with central MDT lesions perform as well as controls for preoperatively learned visual discrimination tasks and in the acquisition of a simple go/no-go odor discrimination task (Slotnick and Kaneko, 1981). However MDT-lesioned rats were impaired when performing odor reversal learning. The authors also noticed in a set of preliminary experiments that three MDT-lesioned rats were not impaired in their acquisition of a visual discrimination reversal set suggesting that the deficits were modality specific. Thus, lesions of the MDT seem to induce a severe deficit in reversal learning and the degree of impairment seems to be related to the extent of the lesion (Slotnick and Kaneko, 1981; Staubli et al., 1987; Lu and Slotnick, 1990; Slotnick and Risser, 1990; McBride and Slotnick, 1997). Importantly, one of the reciprocally connected structures of the MDT—the orbitofrontal cortex—has also been involved in reversal learning (Roesch et al., 2007). The MDT-orbitofrontal cortex network may integrate stimulus-outcome associations to flexibly guide goal-directed behavior.

Rats with thalamic lesions which include the MDT also present deficits in an olfactory continuous delayed nonmatching-to-sample task with no effect on an odor discrimination task. Although, when lesions were more restricted to the MDT, the deficits were minimal in this task (Koger and Mair, 1994; Zhang et al., 1998). The magnitude of the deficits due to MDT lesions is probably related to the extent of the lesions, to the task used, as well as to the difficulty of the task. In fact, Eichenbaum et al. (1980) reported that the major effects of MDT lesions appear when the task complexity is high. Interestingly, the difficulty level of the task may be linked to the attentional demand. An attention deficit may underlie the problems of rodents and humans with damage to the MDT to perform difficult odor discrimination tasks. The role of the MDT in olfactory attention was recently investigated in humans. Plailly et al. (2008) measured attentiondependent network activity using functional magnetic resonance imaging (fMRI). They observed that attending specifically to odor (compared to tone) increases the coupling between piriform cortex to MDT and MDT to orbitofrontal cortex. Corroborating these results, Veldhuizen and Small (2011) observed in humans, using fMRI, activation of the MDT in response to attention to odorants but not tastants. Those two studies support the involvement of the MDT in attention to odors and may be related to the fact that the MDT can be involved in prediction error signaling where the response magnitude of the MDT is significantly higher to unexpected compared to expected odor stimuli (Zelano et al., 2011; Olofsson et al., 2013).

These results obtained with functional imaging are further supported by lesions and neuropsychological studies (Tham et al., 2009, 2011a,b) For example, Tham et al. (2011b) tested whether the MDT was likely involved in top-down directed olfactory attention by using a Target Odor Search Test and showed that patients with damage to the MDT performed poorer verbal-based search than controls. All these studies in humans indicate a possible role of the MDT in olfactory attention processing. However, this idea was debated by Keller (2011) who proposed that the olfactory inputs to the MDT may not be sufficient to support shift of attention toward odors. The question about the contribution of the MDT in olfactory attention still remains open and other studies are required to disentangle it. For example, future studies can assess the impact of lesions of the MDT in rats performing an attention-related task, such as the one described in Ljubojevic et al. (2014) and measure not only their performance relative to controls but also their sampling duration and their latency to reply. Those studies can also assess how and when during the attention-related task the MDT is required by selectively and temporally inhibiting the MDT at different periods of the task using optogenetics.

Lastly, electrophysiological recordings of the MDT in behaving animals may also help to better characterize the temporal contribution of the MDT in olfactory perception and odor-guided behavior. Regarding electrophysiological recordings of the MDT in behaving animals, to the best of our knowledge, there is only one published study recording single-unit activity in the MDT in animals actively engaged in olfactory tasks [in rats: (Kawagoe et al., 2007); nota bene: in Yarita et al. (1980), monkeys were awake though odors were presented passively]. Kawagoe et al. (2007) recorded the MDT activity in an olfactory task requiring animals to discriminate odor cues associated with reward or not. They observed that 10% (13/121 units) of MDT neurons recorded in the central and medial subnuclei responded to the odor cue. Most of these cue-responsive neurons displayed odorant selectivity with a difference of activity between odor cues associated to the same reinforcement category. Combined with this sensory selectivity, the most remarkable effect was that cue-responsive MDT neurons showed strong response preference to cues associated with a reward, and those responses were plastic to extinction and relearning and always related to the reward contingency. The MDT is thus sensitive to stimulus-reward association. Interestingly, the basolateral amygdala and the olfactory tubercle, which project to the medial and central subnuclei, respectively, can also encode the associated outcome of odors, which could contribute to this MDT activity (Schoenbaum et al., 1998; Gadziola et al., 2015). Finally, we recently recorded MDT units in rats performing a two alternative odor discrimination task (Courtiol and Wilson, 2015) and observed that a subset of units were odor selective. Intermingled with this sensory function, we observed that the MDT units displayed activity prior to odor sampling, presumably anticipatory activity, and seemed to encode the choice directiongoal location associated with the odor. This study as well as the one by Kawagoe et al. (2007) emphasize that the MDT may encode both basic sensory as well as complex olfactory functions. This observation is not unique to the MDT as coding both basic sensory and complex olfactory functions has notably been reported in piriform cortex and orbitofrontal cortex (Schoenbaum and Eichenbaum, 1995; Feierstein et al., 2006). With these limited images of MDT function, many questions remain, most notably what are the specific contributions of the MDT in these functions compared to the piriform cortex or the orbitofrontal cortex and what is the MDT adding to olfactory processing?

### Conclusion and Perspectives

Despite the unusual anatomy of the olfactory pathway, an olfactory thalamus can be identified—the MDT receives direct input from various olfactory areas (**Figure 1A**). By virtue of this specific pathway, the thalamic contributions to olfaction are woefully unexplored. Lesion studies in the 1980's followed by more recent humans studies have provided the first evidence of the involvement of the MDT in olfactory processing and suggest a role for the MDT in functions ranging from olfactory perception to attention. While this work provides a glimpse of the place of the MDT in olfaction, many questions, as described above, remain unanswered. For example: Is the MDT role in olfaction similar to the role of primary sensory relays such as the lateral geniculate nucleus? If so, can principles of thalamic function be generalized to all sensory systems? We hypothesize that not to be the case. In olfaction, the functions of ''primary sensory thalamic relay'' including sensory coding, gain control, and state-dependent modulation may be distributed between the olfactory bulb and the piriform cortex (Murakami et al., 2005; Kay and Sherman, 2007). Moreover, given: (1) the anatomical place of the MDT in the olfactory pathway (convergence of many olfactory inputs); (2) the implication of the MDT in higher order functions including olfactory attention; and (3) the MDT being a higher order thalamic relay, the MDT may then be viewed as a higher order olfactory thalamus rather than a primary sensory thalamic relay. Identifying the precise role of MDT in olfactory perception and odor-guided behavior may be an excellent avenue for exploring broader questions of higherorder thalamic function.

### References

### Acknowledgments

This work was supported by grants R03DC014540 from the NIDCD to EC and R01DC003906 from the NIDCD to DAW.

**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 Courtiol and Wilson. 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.

# Efficient population coding of naturalistic whisker motion in the ventro-posterior medial thalamus based on precise spike timing

Michael R. Bale1,2,3 \* † , Robin A. A. Ince2,4† , Greta Santagata2 and Rasmus S. Petersen2 \*

1 School of Life Sciences, University of Sussex, Brighton, UK, 2 Faculty of Life Sciences, University of Manchester, Manchester, UK, 3 Instituto de Neurociencias Alicante UMH-CSIC, Sant Joan d'Alacant, Spain, 4 Institute of Neuroscience and Psychology, University of Glasgow, Glasgow, UK

#### Edited by:

Vincenzo Crunelli, Cardiff University, UK

#### Reviewed by:

Ileana Hanganu-Opatz, Center for Molecular Neurobiology Hamburg (ZMNH), Germany Paloma Tamara Gonzalez-Bellido, University of Cambridge, UK

#### \*Correspondence:

Michael R. Bale, School of Life Sciences, University of Sussex, CRPC Building, Brighton BN1 9QG, UK<EMAIL_ADDRESS>Rasmus S. Petersen, Faculty of Life Sciences, University of Manchester, Stopford Building, Manchester M13 9PT, UK<EMAIL_ADDRESS>
†Joint first authors.

Received: 12 June 2015 Accepted: 03 September 2015 Published: 25 September 2015

#### Citation:

Bale MR, Ince RAA, Santagata G and Petersen RS (2015) Efficient population coding of naturalistic whisker motion in the ventro-posterior medial thalamus based on precise spike timing. Front. Neural Circuits 9:50. doi: 10.3389/fncir.2015.00050 The rodent whisker-associated thalamic nucleus (VPM) contains a somatotopic map where whisker representation is divided into distinct neuronal sub-populations, called "barreloids". Each barreloid projects to its associated cortical barrel column and so forms a gateway for incoming sensory stimuli to the barrel cortex. We aimed to determine how the population of neurons within one barreloid encodes naturalistic whisker motion. In rats, we recorded the extracellular activity of up to nine single neurons within a single barreloid, by implanting silicon probes parallel to the longitudinal axis of the barreloids. We found that play-back of texture-induced whisker motion evoked sparse responses, timed with millisecond precision. At the population level, there was synchronous activity: however, different subsets of neurons were synchronously active at different times. Mutual information between population responses and whisker motion increased near linearly with population size. When normalized to factor out firing rate differences, we found that texture was encoded with greater informational-efficiency than white noise. These results indicate that, within each VPM barreloid, there is a rich and efficient population code for naturalistic whisker motion based on precisely timed, population spike patterns.

#### Keywords: vibrissa, VPM, neural coding, information theory, synchrony, somatosensory

### Introduction

Thalamus is the gateway to the cerebral cortex and the manner in which sensory signals are encoded by thalamic populations fundamentally constrains cortical computation (Sherman and Guillery, 2009). Although it has long been suspected that neuronal activity within modular, neuronal populations is crucial for thalamo-cortical processing, few studies have simultaneously recorded the activity of multiple neurons within such modules in the thalamus (Desbordes et al., 2008; Temereanca et al., 2008; Wang et al., 2010b). The whisker system is ideal for investigating thalamic population coding (Petersen et al., 2009), since there is a well-defined, modular population of neurons devoted to each whisker at each of the major stations of the sensory pathway (Woolsey and Van der Loos, 1970). In the ventroposterior medial nucleus (VPM) of the thalamus, each whisker is primarily represented by a cluster of ∼250 neurons (in rat) known as a ''barreloid'' (van der Loos, 1976).

In the whisker system, mechanoreceptors embedded in the whisker follicle-sinus complex project to the cerebral cortex via parallel, topographic pathways through brainstem and thalamus (reviewed by Diamond and Arabzadeh, 2013). The major (lemniscal) pathway, primarily responsible for encoding whisker-object touch, involves VPM (Deschênes et al., 2005). The core of VPM receives ascending input from the principal nucleus of the brainstem and projects primarily (although not exclusively) to layer IV of primary somatosensory cortex (S1) with collaterals to the thalamic reticular nucleus (NRT). Analogously to other sensory systems, VPM conveys ascending sensory signals to cortex and modulates those signals depending on cortical and behavioral state (Temereanca et al., 2008; Sherman and Guillery, 2009). VPM activity is modulated by input from S1 layer 6, NRT and brainstem neuromodulatory centers. VPM is neurochemically more homogeneous than other thalamic nuclei, such as LGN, in being devoid of GABAergic interneurons and consisting of glutamatergic, S1-projecting relay neurons.

By recording the activity of individual neurons, previous studies have determined that VPM neurons can respond to whisker deflection with high (sub-millisecond) spike timing precision and have identified sensory features to which VPM neurons respond (Waite, 1973; Simons and Carvell, 1989; Pinto et al., 2000; Brecht and Sakmann, 2002; Castro-Alamancos, 2002; Yu et al., 2006; Montemurro et al., 2007a; Petersen et al., 2008; Ego-Stengel et al., 2012). At the population level, it is known that sudden whisker contact (e.g., a ramp-and-hold stimulus) evokes a response where different neurons fire within a few milliseconds of each other—commonly termed ''synchrony''—and that the degree of synchrony varies with stimulus features (Pinto et al., 2000; Temereanca et al., 2008; Wang et al., 2010b). However, it is not known how VPM neurons, at either single neuron or population level, respond to naturalistic whisker motion. Coding principles for naturalistic stimuli can differ substantially from those for artificial stimuli (Mainen and Sejnowski, 1995; de Ruyter van Steveninck et al., 1997). Our aim here was to investigate how the population of neurons in an individual barreloid of the VPM nucleus encodes texture-induced motion of the corresponding whisker and to test whether its population response is synchronous.

To address these issues, we used multi-microelectrode arrays to record simultaneously the activity of multiple single units from within a single VPM barreloid, in response to play-back of texture-induced whisker motion (Arabzadeh et al., 2005; Wolfe et al., 2008; Lottem and Azouz, 2011; Bale et al., 2013). We found that naturalistic whisker motion evoked precisely timed patterns of population spiking, where different groups of neurons fired synchronously at different times (''dynamic synchrony, DS''). Furthermore, mutual information scaled linearly with population size, suggesting a marked lack of redundancy, and responses exhibited a greater informational efficiency for naturalistic stimuli.

## Materials and Methods

### Electrophysiology

All experimental procedures were approved both by the UK Home Office and by The University of Manchester Ethical Review Committee. Adult male Wistar rats (n = 13, 250–350 g) were anesthetized with urethane (1.5 g/kg body weight) and placed into a stereotaxic apparatus. A craniotomy was made over the hemisphere ipsilateral to whisker stimulation 2.0–4.5 mm posterior and 1.5–4.0 mm lateral to bregma. A single-shank, 32-site silicon probe (single row of recording sites, each area 177 µm2 , 50 µm spacing; recording site impedance lowered by Iridium oxide activation, top and bottom sites used to pass current for marking recording site; Neuronexus, Ann Arbor, MI) was inserted into the brain at 45◦ to the midline in the coronal plane (dorso-medial to ventro-lateral), to target the contralateral VPM nucleus. During the experiment, position of the probe within VPM was verified by short-latency responses evoked by manual stimulation of the whiskers. The probe was lowered until whisker-evoked activity ceased, and receptive fields characteristic of the ventroposterior lateral nucleus (VPL) emerged on the deepest recording sites (mean 8.3, SD 0.3 mm from pial surface). Extracellular signals were preamplified, digitized at 24.4 kHz, band pass filtered (300–3000 Hz) and continuously stored to hard disk for offline analysis.

### Whisker Stimulation

All whiskers were trimmed to 5 mm and receptive fields were mapped for each recording site by manually deflecting single whiskers. The whisker that evoked responses at the greatest number of recording sites was selected for stimulation. This whisker was inserted into a pulled pipette tube attached to a piezoelectric actuator (P/N PL127.10; Physik Instrumente) positioned 1 mm from the face. The actuator was shortened to raise resonant frequency and fixed to an aluminium tube via small plastic screws. The dynamic range of the actuator was 0.8 mm (whisker deflection range ∼40◦ ).

We employed two types of dynamic whisker stimulation: lowpass filtered (300 Hz) white noise (Petersen et al., 2008) and a naturalistic texture stimulus (Bale et al., 2013). The texture stimulus was constructed from optical measurements of whisker motion made from awake rats by Wolfe et al. (2008). Wolfe et al. (2008) trained rats to whisk a textured surface (sandpaper: P150, P400, P800 or P1200) and used a CCD array to measure traces of texture-induced whisker motion in the rostro-caudal plane. In order to construct the stimulus corresponding to a given grade of sandpaper, we considered only periods of the CCD traces during which whiskers were in contact with that sandpaper. Such traces were stitched together so that the final position of one trace equalled the first position of the subsequent one until a sequence of 10 s duration was obtained. To minimise potential contribution of head movement, the resulting sequence was highpass filtered (FIR at 1 Hz). Finally, to avoid mechanical resonance from the actuator, the sequence was low-pass filtered (FIR 600 Hz). This procedure was repeated for each of the four grades of sandpaper.

The stimulus protocol consisted of 50–100 trials. Each trial consisted of a fixed 10 s sequence of white noise and 1–4 sequences of 5–10 s repeated texture motion. Stimulus sequences were presented in a randomised order and were separated by 1 s intervals without whisker motion. A schematic for an example session (in which a single texture sequence was presented) is shown in **Figure 2A**.

For consistency with the conditions under which the whisker motion data were registered, all whisker motion stimuli were delivered in the rostro-caudal direction. Accurate reproduction of stimulus position, velocity and acceleration was confirmed by optical testing with a phototransistor circuit (Storchi et al., 2012).

#### Histology

At the end of each experiment, two lesions were made by passing 20 µA for 6 s through both the deepest and most superficial sites on the silicon probe. As detailed in Montemurro et al. (2007a), animals were perfused transcardially with saline and formalin. Brains were removed, immersed in fixative for at least 1 day and then transferred to a phosphate-buffered 30% (w/v) sucrose solution for a further 48 h. Coronal sections of 50 µm were stained with cresyl violet. Recovery of the lesion sites confirmed that all single units reported in this study were located in VPM.

### Spike Sorting

The first step of our data analysis was to isolate single unit activity from the extracellular recordings. We found that spikes emitted by a given single unit were typically recorded from 1–3 recording sites (three recording sites span 100 µm). To exploit this information for better single unit isolation, and also to prevent potential double counting of spikes on adjacent channels, we spike sorted 1–3 channels simultaneously (**Figures 1D–F**). We extracted 1 ms segments around voltage-threshold-crossing times, concatenating segments from 1–3 adjacent recording sites and clustered these data in the space of their principal components using a t-distribution mixture model (**Figure 1E**, Shoham et al., 2003; Bale and Petersen, 2009). Only clusters exhibiting a clear refractory period were accepted. To test for potential double counting of spikes isolated from adjacent blocks of recording sites, we checked for suspicious peaks (width 0.08 ms, corresponding to two sampling intervals) in the cross-correlation function (bin size = 0.04 ms). Any units contributing to such peaks were excluded. An example recording of four simultaneously recorded adjacent channels is shown in **Figure 1D**, of which two sites were sorted together (channels 5 and 6). The spikes extracted during the sorting routine (from **Figure 1F**) are overlaid.

#### Localization of Receptive Fields

During each recording, we identified the whisker that evoked the greatest multi-unit activity (MUA) on each recording site (its principal whisker, PW) to manual deflection of individual whiskers. To assess the multi-unit PW quantitatively, we then recorded responses to piezoelectric deflection of the PW in both caudal and rostral directions with a ramp-and-hold protocol. This was repeated for each of the surrounding whiskers. In offline analysis, we computed the MUA spike count evoked by caudal and rostral deflection at each recording site (as detailed above, time window 20 ms), for each whisker. For each site, we determined the whisker that evoked the greatest MUA response. Provided that this was significantly greater than the spontaneous spike count (time window 20 ms; Wilcoxon signed ranks), this was deemed the PW. By repeating this analysis for each recording site, the PW at each site with significant whisker response was obtained (**Figure 1C**).

The majority of single units (91%) had a receptive field that matched the multi-unit receptive field of the channel to which it was most closely located. The remaining 9% of single units had a receptive field that corresponded to the MU activity of an adjacent whisker.

#### Spike Train Sparseness

Sparse codes are metabolically efficient, in that much information is conveyed with few spikes, and are prominent in the brain (Simoncelli and Olshausen, 2001). We measured temporal sparseness (''lifetime sparseness''), using a slight variant of a standard measure (Rolls and Tovee, 1995):

$$S = 1 - \frac{\langle r\_t \rangle^2}{\langle r\_t^2 \rangle} \tag{1}$$

Here rt is the value of the peri-stimulus time histogram (PSTH) in time bin t and the angled brackets denote the mean over time bins (bin size 10 ms). If rt is constant (non-sparse response), S = 0. In contrast, if rt is zero except for a few bins (highly sparse response), S tends to 1 in the limit of many bins. The original Rolls-Tovee measure, a, is related to the one used here through the relation S = 1−a.

#### Stimulus Sparseness

An important characteristic of sensory stimuli is whether or not they are Gaussian-distributed. This is usefully assessed by measuring ''sparseness'' through the following, standard index (Hyvärinen et al., 2009). Given samples x of a sensory signal:

$$S\_{stim} = -\sqrt{\frac{\pi}{2}} \frac{\langle |\mathbf{x} - \boldsymbol{\mu}| \rangle}{\sigma\_{\mathbf{x}}} \tag{2}$$

Where µ is the mean, σx is the standard deviation of x. Sstim = −1 for Gaussian distributions; S >−1 for distributions with a thinner peak around the mean and fatter tails (''sparse'').

#### Unit Responsiveness

To determine whether a given unit was responsive to the stimulus (white noise or texture), the unit's spontaneous firing rate was measured as the spike count in the 0.5 s interval prior to stimulus onset and its evoked firing rate was measured in each successive 0.5 s interval throughout the course of the stimulus. A unit was classified as ''responsive'' if its firing rate in any of the poststimulus windows was significantly higher than that in the prestimulus window (Wilcoxon signed ranks, p < 0.001; Bonferroni corrected for multiple comparisons). Non-responsive units were not considered in later analysis.

#### Jitter Analysis

Since spike timing precision is a limit on a neuron's capacity to transmit information, it is necessary to quantify it (Petersen et al., 2009). To quantify spiking precision, we measured the trial-to-trial variability in spike timing (''jitter'') for each unit, as previously described (Montemurro et al., 2007a; Bale et al., 2015). First, we divided the evoked response into 0.8 ms time bins and averaged over trials to form a PSTH, which was then smoothed with a Gaussian filter (SD 1.6 ms). Firing episodes corresponded to local peaks in the PSTH. To focus on reliable episodes we selected peaks which satisfied the following two conditions for the jitter calculation: (1) firing rate was at least half of the maximum for that unit and (2) firing rate exceeded a threshold set according to the null hypothesis that the unit fired randomly at the same average rate (p = 0.001). To establish the second condition, we repeatedly simulated random spike trains from a homogeneous Poisson process with the same rate as the timeaveraged firing rate of the considered unit. The threshold was set as the 99.9th percentile of these peak heights. For each peak meeting both criteria, we extracted all spikes fired within ±10 ms of the time of the peak and computed spike times relative to the peak time. These time differences were pooled across both trials and peaks and the resulting distribution fitted to a Gaussian. Jitter was defined as the SD of this Gaussian. A minority of units were only weakly modulated by the whisker stimulus and did not exhibit well-defined firing episodes, resulting in high and unreliable values of calculated jitter. In the jitter plot (**Figure 2E**), we do not plot the points which lie outside the whiskers (defined as median ± 2.7 SD). There were seven such points for white noise and 10 for texture out of a total of N = 65 single units.

### Synchrony Analysis

Synchrony is an interesting property of neuronal population activity, since it can enable reliable transmission of information even if individual synapses are weak (Bruno and Sakmann, 2006). As a first test of synchrony, for each pair of simultaneously recorded units, we computed the cross-correlation function for the stimulus-evoked activity on each trial and then averaged the data across trials (binsize 1 ms). This was done separately for the white noise and texture stimuli. For each pair of units, we calculated the peak value of the trial-mean crosscorrelation function (the ''CCG peak'') as well as its standard error. A unit pair was defined as having a significant crosscorrelation peak if the CCG peak exceeded at least three standard errors.

To test for synchrony in our simultaneously recorded units, we registered the response of each unit in 10 ms time bins, and computed the number of bins in which n = 2, 3, 4 . . . units all fired. As a control, we compared this to the number of bins in which such events would be expected by chance from a statistically independent population of units, with the same single unit spike count distributions. To do this, we shuffled the response of each unit independently across time bins. This preserved the overall firing rate and response diversity, while removing the temporal structure. We then computed the same synchrony measure (number of bins in which different numbers of multiple units fired) for this surrogate data set.

With the simplest (static) type of synchrony, each evoked synchronous population response has the same distribution of population response words (that is, the same set of units fires together). Alternatively, different evoked response events might have different distributions of response words. We term this dynamic synchrony (DS), since events at different times evoke different patterns of synchrony. To quantify the extent to which synchrony is dynamic we used an information theoretic approach. We first isolated synchronous population firing events by identifying time bins (20 ms) where two or more units fired on at least 5% of trials. We quantified the degree of DS as the mutual information between the population response words and the peri-stimulus time of the the response event [Ish information estimator (Montemurro et al., 2007b), Panzeri-Treves bias correction (Panzeri and Treves, 1996)]. To normalize for the potential effect of events of different strengths eliciting the same patterns of synchronous activity but with different probabilities, we consider in the information calculation only trials where at least one unit fired. This is equivalent to removing the response word [0, . . ., 0] from the probability distribution and renormalizing appropriately. If all synchronous population events elicit the same pattern of synchronous activity the information computed in this way should be zero. If it is not, this shows that different synchronous population events evoke different distributions of response patterns—hence the synchrony is dynamic. This can be interpreted as a model comparison between a model where all events elicit population response words from the same distribution, and one in which each event elicits responses from a different distribution—a different pattern of synchrony. To determine the significance of the observed information values we use a permutation testing approach. The null hypothesis for this test is that the observed synchrony is static; the same distribution of response words is evoked by each event (albeit with possibly different pooled firing rate). We therefore performed 1000 permutation calculations in which the identity of the events that elicited each response were randomly shuffled. Mutual information values are affected by the number of trials as well as the size of the spaces considered. Since different recordings had different numbers of trials, cells and detected events it is therefore difficult to compare information values directly. To address this we standardized (z-scored) the observed information for each recording with respect to the corresponding set of permutation values. This gives a measure of the statistical strength of the result in favor of DS, and provides a measure in standardised (z-score) units that can be compared across different experiments.

### Mutual Information Analysis

Mutual information is a powerful, non-linear measure of the correlation between two variables. Applied to neural coding, it quantifies the intuitive notion of the ''information'' that a neuron (or neurons) conveys about a stimulus. A key advantage is that it makes minimal assumptions on the functional relationship between the encoded parameter (e.g., whisker position) and the neuronal response; a second advantage is that it provides a rigorous yardstick for comparing alternative, candidate neural codes. Mutual information can give useful physiological insight concerning whether, for example, precise spike timing increases the capacity of a neuron to transmit sensory messages or whether different neurons in a population convey the same messages (redundancy) or complementary ones. Application of information theoretic methods to the whisker system has been reviewed elsewhere (Petersen et al., 2009; Ince et al., 2010).

The aim of this analysis was to quantify how much total information a given population response code conveys about the dynamic whisker stimulus. To this end, we used the spike trains evoked by the repeated whisker motion sequences. We employed the well-established method of Strong et al. (1998) and conducted the analysis for simultaneously recorded populations of single units. For a population of units i = 1, . . ., N, we computed the response **r** as the N-element binary word (r1, r2, . . ., rN) within the time bin (t, t + 10 ms). Here, ri denotes the (binary) response of unit i where 0 indicates that the unit was inactive (no spikes emitted) in that response bin, and 1 indicates that the unit was active (at least 1 spike emitted); t denotes time with respect to stimulus onset (t = 0). In this approach, these time bins are each considered as the response to a particular stimulus; the stimulus set S is then the set of response bins (implicitly, the small segments of dynamic stimulus preceding each response bin). In this way the explicit consideration of particular stimulus features is replaced by the temporal exploration of a range of inputs from an approximately ergodic process. The method therefore takes into account, without any prior assumptions, all possible stimulus features that might be driving the neural response. The particular number of stimuli available depends on the length of the repeated stimulus segment considered as well as the bin size used. The resulting information value quantifies how discriminable, on average, different response bins (and hence the different preceding stimulus segments) are from each other.

The measurements of r were used to estimate the probability P(r|t) of the neural response r in time bin, t, and the probability P(r), the mean of P(r|t) over all time bins. When a neural response is reliably modulated by a stimulus, it evokes similar responses across repeated trials at a given time, and different responses across different time bins: thus, P(r|t) varies systematically with t. In contrast, when a neuron is insensitive to a stimulus, similar responses are produced for all time bins and P(r|t) is similar to P(r). I(R; S) quantifies how much, on average, P(r|t) differs from P(r):

$$I(R;S) = \left\langle \sum\_{r \in R} P(r|t) \log\_2 \frac{P(r|t)}{P(r)} \right\rangle\_t \tag{3}$$

Here the angled brackets denote an average over all stimulus time bins. I(R; S) quantifies how well, on average, an ideal observer could decide which stimulus time bin was presented from observation of the neural response r on a single trial. I(R; S) has units of bits; one bit of information indicates that, on average, the uncertainty about which stimulus bin was presented is reduced by a factor of two after observation of a single response.

For each experimental session consisting of N units, we computed information values for all possible subpopulations of every size, including single unit information values. For each (sub)-population, the sum of these single unit information values for the population members gives the amount of information that would be conveyed by if the units were transmitting information independently (Averbeck et al., 2006), denoted independent. We measured the synergy/redundancy as the difference between the information conveyed by the population and the sum of the single unit information values (Schneidman et al., 2003; Latham and Nirenberg, 2005):

$$\text{Sym}(R\_1, R\_2; \mathcal{S}) = I(R\_1, R\_2; \mathcal{S}) - I(R\_1; \mathcal{S}) - I(R\_2; \mathcal{S}) \tag{4}$$

Positive values of Syn(R1, R2; S) indicate synergy; the two response variables carry more information about the stimulus together than they do alone. Negative values of Syn(R1, R2; S) indicate redundancy; the two variables together carry less information then the sum of their individual contributions. To test whether a low synergy/redundancy value for a given population might reflect cancelation between strongly redundant and strongly synergistic unit pairs, we conducted the following analysis. For a given population, we computed the synergy/redundancy of all of its constituent unit pairs. Positive (synergistic) and negative (redundant) values were summed separately, resulting in a net pairwise synergy value and a net pairwise redundancy value. The maximum absolute value of these two was taken an indicator of the maximum possible pairwise interaction effect for the population. We then normalized this value with respect to the mutual information conveyed by the population.

To give a more intuitive measure for the mutual information conveyed by a population, we converted bits to bits/s by dividing by the bin width. We also separately normalized for differences in firing rate, by calculating the mutual information per spike (bits/spike), both for single cells and populations.

The probability distributions that appear in equation (3) were computed from a finite number of trials and thus subject to sampling error, which typically inflates estimates of mutual information (Panzeri et al., 2007). Bias correction was performed using the shuffled information estimator, Ish(R; S) together with quadratic extrapolation (Montemurro et al., 2007b; Panzeri et al., 2007). Mutual information values were computed using the PyEntropy package (Ince et al., 2009). The length of the repeated stimulus blocks was either 5 s or 10 s, resulting in 500 or 1000 response bins for each information calculation. The crucial parameter for determining the effectiveness of bias correction is the ratio of the number of trials per stimulus to the number of different response symbols (Ince et al., 2010). This number of trials per stimulus is given by the number of repeated presentations of the continuous stimulus and does not change across different stimulus block lengths or response bin sizes. The number of response symbols (2 to the power of the population size) is also fixed, so the difference in stimulus block lengths does not affect the bias of the information measure. With the available data, the number of trials per stimulus was 50–200 (median 200); the ratio of trials/stimulus to number of response symbols (for the largest population size recorded in each experiment) was 0.2–100 (median 18.75). Simulations have shown (Ince et al., 2009) that with the bias correction employed here, a ratio of ∼0.25 can be adequately bias corrected: the data here fall within that range. To further test the accuracy of our bias correction, we repeated the analysis with a surrogate data set in which, for each trial, the population responses were shuffled independently across time bins (so that the true information is zero). The resulting mean information values for each population size were expressed as a percentage of the corresponding unshuffled value. In the most demanding case—population size 9—for white noise stimulation, the median was 0.79%; for texture stimulation, the median was 0.85%. This indicates that bias was effectively corrected.

### Results

### Recording the Population Response of a VPM Barreloid to Whisker Motion

Our primary aim was to investigate how the population of neurons within a single VPM barreloid collectively encodes naturalistic whisker motion. To address this, it was essential to deliver multiple, identical repeats of controlled, whisker motion sequences. To this end, we implanted 32-channel silicon probes into the VPM of urethanised rats, and recorded simultaneously the responses of multiple single units to deflection of the whiskers. To concentrate recording sites as far as possible within the same barreloid, we implanted probes at a 45◦ angle (dorso-medial to ventro-lateral), approximately parallel to the long axis of the barreloids (**Figures 1A,B**). To assess microelectrode placement, we identified, for each recording site, the whisker that evoked the greatest multi-unit response to ramp-and-hold deflection (**Figure 1C**; see ''Materials and Methods'' Section)—the ''PW''. In the example of **Figure 1**, the PW was E2 for sites 4–11, E1 for site 13 and D2 for sites 14–16. Overall (N = 16 recordings), up to 11 recording sites shared the same principal whisker (mean 7.4, SD 2.4 range 4–11), indicating that they recorded neural activity from the same barreloid. Since the recording sites were spaced at 50 µm intervals along the probe shank, 11 sites spanned 500 µm. This length is consistent with anatomical estimates of barreloid length in adult rats (Haidarliu and Ahissar, 2001).

A piezoelectric actuator was used to ''play back'' sequences of dynamic whisker motion (Bale et al., 2013; see ''Materials and Methods'' Section). A naturalistic stimulus was generated from optically registered whisker sweeps of rats actively whisking sandpaper (Wolfe et al., 2008). To benchmark the coding efficiency of this stimulus, we also used low-pass filtered white noise (hereafter abbreviated to ''white noise''; **Figure 2A**), which was normalized to have the same standard deviation. Both stimuli were applied to individual whiskers (N = 16 recordings), and we recorded neuronal responses from all probe sites simultaneously.

### Response of Single VPM Units to Naturalistic Whisker Motion are Temporally Precise and Sparse

**Figure 2** shows PSTHs, evoked by white noise (**Figure 2A**) and texture (**Figure 2B**) applied to whisker E2, for a set of 9 responsive single units (see ''Materials and Methods'' Section) recorded simultaneously from barreloid E2. The single unit responses to the two stimuli were qualitatively similar. Both stimuli evoked responses that consisted of temporally isolated firing episodes. We quantified the spike timing precision of these firing episodes by estimating the trial to trial ''jitter'' in spike time (**Figure 2E**; see ''Materials and Methods'' Section). Jitter for both texture and white noise was sub-millisecond (median 0.48 ms, IQR 0.48–1.2 ms and median 0.36 ms, IQR 0.36–0.88 ms respectively). The main differences between the stimuli were that white noise evoked a higher rate of firing events (**Figure 2C**; median 3.0 spikes/s compared to 1.4 spikes/s; signed rank test, p < 1e−7 ), while responses to texture exhibited higher temporal sparseness (**Figure 2D**; median 0.80 compared to 0.61; signed rank test, p < 1e−9 ). To investigate the source of these differences we compared the stimuli. Compared to white noise, the texture stimulus both had a markedly non-Gaussian ''sparse'' distribution (**Figure 2F**) and contained less relative power at high frequencies (**Figure 2G**). Using a standard index which is −1 for a Gaussian and >−1 for sparse distributions (Material and Methods, Equation 2), the white noise had an index value of −1 and the texture −0.85.

### VPM Population Responses Exhibit Dynamic Synchrony

Our finding that texture evokes precisely timed spikes raises the question of the nature of the response at the population level. Motivated by previous studies of the VPM response to periodic whisker stimulation (Bruno and Simons, 2002; Temereanca et al., 2008; Wang et al., 2010b), we investigated whether the population response manifested synchrony, conventionally defined as coincident firing on a time-scale of ∼10 ms (Wang et al., 2010a). Inspection of the recordings revealed qualitative evidence for synchrony. **Figure 3A** shows a detailed view of the responses of 3 simultaneously recorded single units (a subset of those shown in **Figure 2**) evoked by the texture stimulus. Within this 2.5 s excerpt, there were several times at which the stimulus evoked a marked increase in firing rate in 2 or more units (shaded regions).

To quantify the synchrony, we first computed crosscorrelograms (CCGs) for the texture-evoked activity (**Figure 3B**, for units with PSTHs shown in **Figure 3A**). Eighty four neuron pairs (out of a total of 145) exhibited statistically significant CCG peaks (see ''Materials and Methods'' Section). For the significant peaks, the CCG peak amplitude was 0.03 ± 0.11

(color coded). (B) Cross-correlograms for each pair-wise combination of the units of panel (A). (C) Log frequency with which a given number of units fired simultaneously (10 ms bin). Frequency is normalized by the single unit firing rate. Solid line shows results for simultaneously recorded data, dashed line shows control under null hypothesis that each unit fired independently at random (see "Materials and Methods" Section).

coincidences/spike (median ± SD), the CCG (absolute) peak lag was 2.5 ± 5.8 ms and the CCG peak width (full width at half maximum, FWHM) was 4.6 ± 3.6 ms. Results for white noise were similar: 71 pairs exhibited significant CCG peaks; for these peaks, CCG peak amplitude was 0.08 ± 0.27; lag was 2.0 ± 5.6 ms and CCG peak FWHM was 3.3 ± 5.1 ms.

To extend the analysis of synchrony beyond the pairwise CCGs to neuronal sets of arbitrary size, we used the following procedure. We registered the response of each unit in 10 ms time bins, and computed the number of bins in which n = 2, 3, 4 . . . units all fired on the same trial. As a control, we compared this to the number of such events expected from random firing (see ''Materials and Methods'' Section). We found, for the textureevoked response of the nine simultaneously recorded units of **Figure 2A**, that synchronous firing of three and four units occurred more often (14.6 and 66.8 times respectively) than expected from random firing (**Figure 3C**). Synchronous firing of 5 units occurred in the experimental data but never in the control. On average, across all sets of simultaneous recordings, synchrony occurred more frequently than expected by random firing.

These results indicate that both white noise and texture whisker motion induced synchronous firing within the associated VPM barreloid. However, closer inspection of the data revealed that the precise pattern of synchrony changed over time: different constellations of neurons fired at different times during the stimulus. The shaded regions in **Figure 3A** illustrate four response episodes where the pooled firing rate of the population was similar but the mean population spike patterns were, in 3 cases, distinct. For two of the peaks all three units exhibited high firing rate at almost the same time (red shading), but at other times different combinations of units responded: units 1 and 2 but not 3 (blue shading) or units 1 and 3 but not 2 (green shading).

With the simplest (static) type of synchrony, each evoked synchronous population response has the same distribution of population response words (that is, the same set of units fires together). Alternatively, different evoked response events might have different distributions of response words. Since events at different times evoke different patterns of synchrony, we refer to this type of synchrony as dynamic (dynamic synchrony, DS). DS implies that there is a repertoire of different synchronous response patterns and that there is a systematic relationship between a particular stimulus event and the evoked response pattern. DS can be quantified by measuring the strength of this association. Specifically, we quantified DS with a normalized mutual information value (see ''Materials and Methods'' Section). This value can be interpreted as a z-score: under the null hypothesis of static synchrony, it should take values drawn from a standard normal distribution. For all seven recordings in which synchronous events were reliably identified (3–9 simultaneously recorded units; 3–121 synchronous events) the DS measure was highly significant (p < 0.001; permutation test, N = 1000). We found that responses to texture stimuli resulted in a DS measure of 35.8 ± 21.1 (standardised units) while for those to white noise, the values were 26.8 ± 23.1; but that this difference was not significant (Wilcoxon signed rank test, p > 0.3).

### VPM Populations Encode Information Independently

To guide our investigation of the population response, we examined how the information conveyed by the population scaled with population size. Such an analysis can reveal whether the representation is, overall, redundant (sub-linear scaling), independent (linear scaling) or synergistic (super-linear scaling). We estimated the mutual information that the population spike pattern (the simultaneous response of all neurons in a given time bin) conveyed about white noise and texture (see ''Materials and Methods'' Section). In all cases, we found mutual information to scale approximately linearly with population size. We illustrate this with the three experiments with highest numbers of simultaneously recorded units (**Figure 4**). The similarity between the mutual information available in the population and that in the sum of single unit information values, across a range of population sizes (**Figures 4A,B**, independent solid lines, pattern dotted lines), implies a striking lack of redundancy.

In principle, zero redundancy for a population could mask significant redundancy/synergy in subsets of the population that cancel out. To test this, we first calculated pairwise synergies and redundancies directly (see ''Materials and Methods'' Section). For each recording with two or more simultaneous units, we calculated the mean absolute synergy/redundancy value (Equation 4) of each unit pair, as a percentage of the mutual information conveyed by that pair. Across recordings (N = 11), the mean was 2.1 ± 1.7% for white noise and 1.8 ± 1.0% for textures. This shows that net pairwise interaction effects were small. Second, for recordings in which at least five units were recorded (N = 3; those shown in **Figure 4A**), we considered all possible subpopulations of size 5. For each such sub-population, we computed an index which expresses the maximum possible pairwise synergistic/redundant effect, as a percentage of the mutual information conveyed by the sub-population (see ''Materials and Methods'' Section). This resulted in a maximum pairwise synergy/redundancy cancellation of 2.0 ± 1.5% (N = 258 populations, mean ± SD). While this analysis does not exclude the possibility of potential higher order effects, it does show that the near-zero redundancy observed at the population level was not a consequence of cancelation between synergy and redundancy in pairs of units.

FIGURE 4 | Quantitative comparison of alternative population codes. (A) Mutual information conveyed about the white noise whisker stimulus by the VPM population (bin size 10 ms) for the three experiments with the largest number of simultaneously recorded units. "Population" denotes the information in the synchronous population response; "Independent" denotes the linear sum of the information conveyed by each individual neuron within the considered population. Each point is an average over all possible sub-populations of a given size (error bars show SEM over populations). (B) Corresponding results for responses to sandpaper texture stimulus. (C) Information conveyed by the different codes under white noise (blue) and texture (red) stimulation for all subpopulations obtained over all experimental sessions. Information values are normalized by the total firing rate of each considered sub-population to obtain units of bits/spike. Each point is the mean over all sub-populations of that size; error bars show SEM. (D) Scatter plot of raw information values for single units, pairs of units and populations of size 5.

### VPM Populations Encode Naturalistic Texture Motion Efficiently than White Noise Motion

We repeated the information calculations for all recordings (N = 16; **Figure 4D**). To factor out the effect of different firing rates, we normalized the mutual information values to units of bits/spike and then averaged the results over all simultaneous recordings of a given size (**Figure 4C**). The normalization to units of bits/spike revealed that the mutual information conveyed per spike was higher for the texture stimulus than for the white noise stimulus. At the single unit level (N = 65), white noise stimulation yielded mutual information values of 1.33 ± 0.67 [1.29] bits/spike, while texture stimulation yielded 1.94 ± 1.24 [1.64] bits/spike (mean ± SD [median]): this difference was significant (signed ranks test; p < 0.00001). For populations of size 6 (N = 169), the population code conveyed 1.26 ± 0.15 [1.26] bits/spike with white noise and 1.95 ± 0.26 [1.98] bits/spike with texture (signed ranks test; p < 1e−10). This difference might mean either that the population transmits more information (in a given period of time) about texture than about white noise or, alternatively, that it conveys the same information but using fewer spikes. The scatter plots in **Figure 4D** show that the mutual information values conveyed by the response in each 10 ms bin were usually similar for the two types of stimuli (although a proportion of population responses were more informative about white noise). Therefore, the lower firing rate evoked by texture (**Figure 2C**) accounts for the finding that texture conveys more bits/spike. The population response conveys about the same amount of information about texture motion as white noise, but does so with fewer spikes. This suggests that the subcortical pathway may be adapted to whisker motion with the statistical characteristics of natural object interactions.

### Effect of Noise Correlations on the Population Code

Dynamic synchrony might reflect either: (1) properties of individual units or (2) network properties. To discriminate between these possibilities, we calculated the effect of noise correlations on the population code, using information theory (Schneidman et al., 2003; Latham and Nirenberg, 2005; Chicharro, 2014). Noise correlations are correlations between the responses of simultaneously recorded single units that are not related to, or driven by, the external stimulus; arising instead from intrinsic network activity. The effect of noise correlations can be removed by creating a surrogate data set in which the responses of each unit to a given stimulus are shuffled across trials, with respect to the other members of the population. This surrogate data set is then a model for what the responses would look like if the cells were firing independently. We generated such surrogate data and repeated both the DS analysis and the mutual information based analysis. We found, for both types of stimulus, that shuffling made little difference to the DS measure: the ratio of the shuffled to the unshuffled DS measure was 1.02 [0.82 1.19] (median [min max] across recordings) for white noise and 0.99 [0.90 1.13] for texture. To assess the impact of noise correlations on mutual information, we computed the difference between the unshuffled and shuffled information, normalized by the unshuffled information. Again, for both types of stimulus, the impact of noise correlations was small: 4.1% [2.4 6.7] for white noise, 2.9% [1.8 4.8] for textures (median [min max] across population sizes). These results indicate that the results reported above reflect the coding properties of individual units rather than network interactions.

### Robustness of Results to Differential Parameters of the Mutual Information Estimation

It is necessary to test whether estimates of mutual information are accurate (see ''Materials and Methods'' Section). To assess the robustness of the main mutual information results to the various parameters of the information theoretic analyses, we performed the following controls. First, we investigated the effect of the length of stimulus. In the full data set, stimulus lengths ranged from 5–40 s, corresponding to 500–4000 10 ms bins. To test whether our results were affected by these differences, we repeated the population mutual information analysis using subsections of the stimulus sequence (1.25–5 s) common to all recordings (**Figure 5A**). In all cases, the texture stimulus evoked more informative responses than white noise (in bits/spike), there was very little redundancy and the values were quantitatively very close to the full results presented in **Figure 4C**. Information values for the shortest texture stimuli (1.25s) were lower, but we expect this is a consequence of the sparse structure of the stimulus. Second, we tested the effect of bin size used for the information analysis (range 10–80 ms; **Figures 5B,C**). Again the results were robust: in all cases, texture evoked higher bits/spike than white noise and redundancy was very low. Information values were lower for larger bin sizes; this indicates the temporal precision of the population response is (as expected from the jitter and CCG results) at least 10 ms.

### Discussion

The aim of this study was to investigate how a modular, population of neurons corresponding to one whisker processes complex, naturalistic whisker motion. We focussed on the VPM thalamus, where signals from a given whisker are primarily encoded by a population of ∼250 projection neurons within the corresponding barreloid (van der Loos, 1976; Land and Simons, 1985; Sugitani et al., 1990; Haidarliu et al., 2008). We were able to record from up to nine single units within one barreloid simultaneously. In contrast to what has previously been reported with ramp-and-hold stimuli, we found that naturalistic whisker motion evoked a dynamic sequence of population spiking, where different groups of neurons fired synchronously at different times (''dynamic synchrony'').

Previous work on the rat/mouse VPM recorded the activity of 1–2 neurons within one barreloid (Waite, 1973; Simons and Carvell, 1989; Pinto et al., 2000; Brecht and Sakmann, 2002; Yu et al., 2006; Montemurro et al., 2007a; Petersen et al., 2008; Temereanca et al., 2008; Bale and Petersen, 2009; Wang et al., 2010b; Scaglione et al., 2011; Poulet et al., 2012). Here, by using multi-site silicon probes and a novel insertion angle,

up to 11 recording sites could be located within a single barreloid. Eleven sites spanned 500 microns, which is consistent with anatomical measurements of barreloid length in adult rats (Haidarliu et al., 2008). This technique enabled us to isolate the activity of up to nine single units within the same barreloid.

''Naturalistic'' stimuli, that reflect the complex, dynamic stimuli that animals experience during natural behavior, can evoke responses that differ substantially in reliability and spike timing (Mainen and Sejnowski, 1995; de Ruyter van Steveninck et al., 1997). Although the naturalistic stimulus paradigm is wellestablished in vision (Simoncelli and Olshausen, 2001), it is only more recently that it has become feasible in the whisker system (Arabzadeh et al., 2005; Lottem and Azouz, 2011; Bale et al., 2013). Our aim here was to study the population response of a VPM barreloid to naturalistic whisker motion. To this end, it was essential to reproduce identical sequences of naturalistic whisker motion on multiple trials. Our approach was to use an anesthetised preparation and to play back sequences of textureinduced whisker motion recorded optically from behaving rats (Wolfe et al., 2008).

Previous work on VPM has shown that rapid whisker motion, such as occurs frequently during a white noise sequence or at the start of a ramp-and-hold stimulus, evokes spikes whose timing is reliable with sub-millisecond precision (Montemurro et al., 2007a). Under some circumstances (periodic whisker deflection), it is known that these spikes occur coincidently across neurons on a time-scale of ∼10 ms (''synchrony''; Bruno and Sakmann, 2006; Temereanca et al., 2008; Wang et al., 2010b). We found that naturalistic whisker motion evoked a complex population response where different subsets of neurons are co-active at different times (**Figure 3A**). In other words, the synchrony was dynamic. We argue that this difference is due to the more complex nature of the stimuli we used. The onset of a ramp-and-hold whisker deflection provokes a correlated increase in position, velocity, acceleration and all higher order temporal derivatives of whisker position (Petersen et al., 2008). Thus, although VPM neurons are tuned to diverse kinetic features (Pinto et al., 2000; Petersen et al., 2008), such a stimulus will tend to evoke a synchronous response from a substantial proportion of the neuron population. Conversely, with a dynamic naturalistic stimulus, the position, velocity and acceleration are decoupled and neurons tuned to different kinetic features tend to fire at different times. Consistent with this interpretation, our quantitative measure of DS was not strongly affected by removing noise correlations from the responses. Moreover, the degree of synchrony exhibited by VPM neuron pairs is known to depend on both deflection velocity and stimulation frequency (Temereanca et al., 2008). Adaptation mechanisms may contribute further dynamism to the response (Wang et al., 2010b). VPM response heterogeneity has also been reported both in the electrical whisking paradigm (Yu et al., 2006) and in the behaving mouse (Gutnisky et al., 2013). These observations are consistent with a model of VPM barreloid function whereby different neurons independently encode different aspects of the on-going whisker motion. Consistent with this, we found that the mutual information conveyed by the population scaled linearly with population size and that this effect was robust to varying the parameters of the analysis. We found no evidence that this lack of redundancy reflected cancelation between strongly synergistic and strongly redundant pairwise sub-populations.

Individual spikes conveyed more bits/spike for the naturalistic texture stimulus compared to white noise (**Figure 4C**), suggesting that the sensory pathway to thalamus may be optimized towards representing naturalistic stimuli, as has also been suggested in the visual and auditory systems (Dan et al., 1996; Simoncelli and Olshausen, 2001; Lewicki, 2002). Similarly to these other systems, the naturalistic texture stimulus had a notably non-gaussian ''sparse'' distribution. Compared to white noise, the texture stimulus had relatively greater power at low frequencies. There are likely to be two factors underlying the greater sparseness and informational efficiency of the texture responses: (1) the greater stimulus sparseness and (2) the fact that VPM neurons tend to be velocity-sensitive (Pinto et al., 2000; Petersen et al., 2008) and therefore sensitive to higher stimulus frequency components. Theories of efficient coding (redundancy reduction, predictive coding) predict that, under high signal to noise conditions, there should be little redundancy in the neural representation of a sensory signal (Barlow, 1961; Srinivasan et al., 1982; Atick and Redlich, 1990). There has been extensive investigation of efficient coding in vision (Atick and Redlich, 1990; Simoncelli and Olshausen, 2001; Vinje and Gallant, 2002; Sharpee et al., 2006), a few studies in audition (Smith and Lewicki, 2006; Ming and Holt, 2009) but in somatosensation, the only previous work has been on tactile robots (Hafner et al., 2003; Evans, 2013). As far as we are aware, the current study is the first direct evidence for efficient coding in somatosensation.

Synchrony is potentially important, since it provides a possible explanation for how the thalamocortical pathway might achieve reliable transmission, despite the fact that individual, afferent synapses are typically weak (Bruno and Sakmann, 2006). When a population of thalamic neurons spikes synchronously (on the time-scale of membrane time constant), multiple postsynaptic potentials can summate in a thalamorecipient neuron and trigger postsynaptic spiking. Such an integration mechanism accounts for how stimuli such as ramp-andhold whisker deflections might trigger a reliable response in thalamorecipient, cortical neurons, and is a potentially effective method of decoding the occurrence and velocity/direction of such stimuli. However, integration is sensitive only to the number of presynaptic spikes and is blind to their origin. It will lose information from a DS code where different patterns of spikes are elicited by different stimulus features. An intriguing possibility is that, through tuning of synaptic strengths (Feldman, 2009) and/or nonlinear dendritic processes (London and Häusser, 2005), cortical neurons are able to perform decoding operations more sophisticated than integration and thereby to discriminate amongst a broad class of sensory signals.

### References

Arabzadeh, E., Zorzin, E., and Diamond, M. E. (2005). Neuronal encoding of texture in the whisker sensory pathway. PLoS Biol. 3:e17. doi: 10.1371/journal. pbio.0030017

Our findings raise issues for further work. One limitation of the current analysis relates to the information calculations. Mutual information as computed here is a powerful statistic for quantifying, in a general way, the reliability of the population response to a dynamic stimulus. However, it would also be interesting to understand what particular kinematic features of the stimulus drive the population response. Second, although recording from populations of size 9 from a single barreloid is the largest-scale recording yet reported from VPM, some network interaction effects only become appreciable for large population sizes (Roudi et al., 2009), and it would be interesting to achieve still higher population sizes. Finally, it is an important challenge to investigate population coding in the awake, behaving animal. However, since multiple trials of identical whisker motion sequences cannot be delivered under these conditions, the present analysis approaches are inapplicable and new methods will be required.

In sum, our findings suggest that the basic building block of the whisker-related thalamus—the barreloid—encodes naturalistic sensory information in a remarkably efficient manner. Barreloid ensembles exhibit temporally sparse responses that encode rich information about whisker motion through a DS code. The combination of sub-millisecond spike timing precision and diverse tuning leads to a high capacity thalamic population code, which potentially conveys a rich signal about multiple features of whisker motion to cortex. A potentially significant implication of our findings is that it would be advantageous for cortical circuits to act not only as detectors of coincident thalamic activity but also to decode patterns of thalamic population activity.

### Author Contributions

MRB, RAAI and RSP designed the experiment. MRB and GS performed neurophysiology. MRB, RAAI and RSP performed data analysis. MRB, RAAI and RSP wrote the manuscript. MRB and RAAI contributed equally and share first authorship of this manuscript.

### Acknowledgments

We thank M. Montemurro and M. Evans for valuable discussions, D. Feldman and J. Wolfe for sharing their whisker motion data and for valuable discussions. This work was supported by the CARMEN e-science project (EPSRC grant EP/E002331/1), BBSRC BB/G020094/1 and BB/L007282/1, MRC MR/L01064X/1, Wellcome Trust 097820/Z/11/B, Spanish Ministry of Science and Innovation (BFU2011-23049, co-funded by the European Regional Developmental Fund) and the Lord Alliance Foundation.

Atick, J. J., and Redlich, A. N. (1990). Towards a theory of early visual processing. Neural Comput. 2, 308–320. doi: 10.1162/neco.19<IP_ADDRESS>8

Averbeck, B. B., Latham, P. E., and Pouget, A. (2006). Neural correlations, population coding and computation. Nat. Rev. Neurosci. 7, 358–367. doi: 10. 1038/nrn1888

**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 Bale, Ince, Santagata and Petersen. 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.

# A computational relationship between thalamic sensory neural responses and contrast perception

Yaoguang Jiang1 , Gopathy Purushothaman2 and Vivien A. Casagrande1,2,3 \*

1 Department of Psychology, Vanderbilt University, Nashville, TN, USA, 2 Department of Cell and Developmental Biology, Vanderbilt University, Nashville, TN, USA, 3 Department of Ophthalmology and Visual Sciences, Vanderbilt University, Nashville, TN, USA

Uncovering the relationship between sensory neural responses and perceptual decisions remains a fundamental problem in neuroscience. Decades of experimental and modeling work in the sensory cortex have demonstrated that a perceptual decision pool is usually composed of tens to hundreds of neurons, the responses of which are significantly correlated not only with each other, but also with the behavioral choices of an animal. Few studies, however, have measured neural activity in the sensory thalamus of awake, behaving animals. Therefore, it remains unclear how many thalamic neurons are recruited and how the information from these neurons is pooled at subsequent cortical stages to form a perceptual decision. In a previous study we measured neural activity in the macaque lateral geniculate nucleus (LGN) during a two alternative forced choice (2AFC) contrast detection task, and found that single LGN neurons were significantly correlated with the monkeys' behavioral choices, despite their relatively poor contrast sensitivity and a lack of overall interneuronal correlations. We have now computationally tested a number of specific hypotheses relating these measured LGN neural responses to the contrast detection behavior of the animals. We modeled the perceptual decisions with different numbers of neurons and using a variety of pooling/readout strategies, and found that the most successful model consisted of about 50–200 LGN neurons, with individual neurons weighted differentially according to their signal-to-noise ratios (quantified as d-primes). These results supported the hypothesis that in contrast detection the perceptual decision pool consists of multiple thalamic neurons, and that the response fluctuations in these neurons can influence contrast perception, with the more sensitive thalamic neurons likely to exert a greater influence.

Keywords: lateral geniculate nucleus (LGN), perception, contrast, neural model, choice probability

### Introduction

From smelling a flower to recognizing the face of a loved one, every perceptual task we face, simple or complex, involves a number of neurons in a wide range of brain areas. Of essential interest to neuroscientists is the number of sensory neurons needed to sustain a perception, and the way these neurons are decoded at a later stage to form various decisions. Theoretically, every perceptual task can be accomplished by engaging only the few sensory neurons that are the most sensitive for that task (i.e., the lower envelope principle; Barlow, 1995; Parker and Newsome, 1998).

### Edited by:

W. Martin Usrey, University of California, Davis, USA

#### Reviewed by:

Donald A. Wilson, New York University School of Medicine, USA Henry Joseph Alitto, University of California, Berkeley, USA

#### \*Correspondence:

Vivien A. Casagrande, Department of Cell and Developmental Biology, Vanderbilt University, PMB407935 RM U-3218 MRB 3, 465 21st Avenue South, Nashville, TN 37240-7935, USA<EMAIL_ADDRESS>
Received: 28 July 2015 Accepted: 14 September 2015 Published: 08 October 2015

#### Citation:

Jiang Y, Purushothaman G and Casagrande VA (2015) A computational relationship between thalamic sensory neural responses and contrast perception. Front. Neural Circuits 9:54. doi: 10.3389/fncir.2015.00054 In reality, however, a variety of factors such as the response variances of single neurons and the positive noise correlations between pairs of neurons constrain the pool size, requiring at least 10–1000 sensory neurons in an average sized decision pool (Shadlen et al., 1996; Cook and Maunsell, 2002; Purushothaman and Bradley, 2005; Cohen and Newsome, 2009; Liu et al., 2013). In sensory cortex, such perceptual decision pools have two prominent features. First, stimulus-independent, random fluctuations of sensory neural responses are known to covary with the perceptual decisions of the animal. The strength of this covariation is quantified as ''choice probability'' (Britten et al., 1996). Weak but significantly above chance choice probabilities have been observed in a number of sensory cortical areas (Britten et al., 1996; Dodd et al., 2001; Cook and Maunsell, 2002; Grunewald et al., 2002; Uka and Deangelis, 2004; Liu and Newsome, 2005; Purushothaman and Bradley, 2005; Uka et al., 2005; Nienborg and Cumming, 2006; Palmer et al., 2007). Second, cortical sensory neurons are also correlated with each other in their random response fluctuations (Averbeck et al., 2006; Cohen and Kohn, 2011). This correlation, known as the interneuronal noise correlation, is likely to reflect the shared feedforward, feedback, or lateral connections between neurons (Zohary et al., 1994; Shadlen and Newsome, 1998; Bair et al., 2001; Reich et al., 2001; Cohen and Maunsell, 2009). Previous modeling work has revealed that interneuronal correlations can have a profound influence on the choice probability structure of the decision pool (Shadlen et al., 1996; Cohen and Newsome, 2009; Nienborg and Cumming, 2010; Haefner et al., 2013).

Such interneuronal correlation or choice probability measurements, however, are rarely made in subcortical structures (but see Liu et al., 2013). In the mammalian visual system, the retina sends direct input to the lateral geniculate nucleus (LGN) of the thalamus which, in turn, relays this information to the visual cortex. Recently, we reported the first study in which LGN neural responses were examined in detail while the animals were required to make perceptual decisions using the information available within the receptive fields of those LGN neurons (Jiang et al., 2015). In a two alternative forced choice (2AFC) contrast detection task, we found that the majority of single LGN parvocellular (P) and magnocellular (M) neurons were not as sensitive as the monkeys. Importantly, the covariation between neural responses and perceptual decisions, measured as choice probability, was significant for both P and M neurons, even though the average interneuronal correlation between LGN neuron pairs was not different from zero. Additionally, both neural sensitivity and choice probability evolved throughout the stimulus presentation time, with M neurons exhibiting faster and more transient response profiles than P neurons (Jiang et al., 2015).

Taking advantage of this previously characterized dataset and using a computational approach, we investigated in this study how single LGN neurons contribute to our perception of contrast. We built a series of models to explore the interaction between the size of the decision pool, the duration of integration time, and the pooling/readout strategy of the neural system. Because previous experimental and computational work has suggested a positive relationship between neural sensitivity, choice probability, and readout weight (for example see Britten et al., 1996; Shadlen et al., 1996; Purushothaman and Bradley, 2005; Haefner et al., 2013; Liu et al., 2013), we examined not only the standard uniform readout model but also several alternative weighted readout schemes in which individual neurons were assigned different weights based on their sensitivities. We accepted or rejected these models based on their ability to account for the behavioral performance of the monkeys as well as the measured choice probability values for LGN neurons (see Jiang et al., 2015). Aspects of the modeling data presented here have been published in abstract form (Jiang et al., 2012, 2013).

### Materials and Methods

All the experimental procedures regarding surgical preparation, animal training, stimulus presentation, and physiological recordings have been described in detail in previous publications (Jiang et al., 2013, 2015), and are therefore only briefly repeated here when relevant.

### Subjects

Two macaque monkeys (monkey 1: Macaca radiata, male, 7 kg, 10 years old; monkey 2: Macaca mulatta, male, 8 kg, 12 years old) served as subjects. The monkeys were treated and cared for in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and the guidelines of Vanderbilt University Animal Care and Use Committee under an approved protocol. The monkeys underwent sterile procedures for the implantation of head posts and recording chambers. The chambers were centered over the right LGN of monkey 1 (AP = 7, ML = 12.5) and the left LGN of monkey 2 (AP = 7, ML = 12).

### Visual Stimulus Presentations and Behavioral Tasks

The monkeys were first trained to fixate on a central fixation spot for an extended period of time. Next, the monkeys were trained to perform a two-alternative forced choice (2AFC) contrast detection task, in which a contrast stimulus was presented either at the receptive field location of the cell being recorded, or at a symmetrical location in the opposite visual hemi-field, for a fixed duration (200 ms). The monkeys saccaded to one of the two target locations to indicate the side on which the stimulus was presented. The stimulus diameter was always the sum of the classical receptive field diameter (center and surround) plus the fixation window diameter. During each recording session, stimuli of 5 or 9 different contrast levels (including 0% contrast, or blank trials, where no physical stimulus was presented) were presented at each location. Different contrast levels and presentation locations (i.e., left or right) were randomly mixed, with equal probabilities of left or right appearance and higher proportions of low to medium contrast trials to ensure accurate estimations of the psychophysical threshold.

### Psychometric Functions

The proportion of correct responses from the monkeys was plotted for each contrast, and a Weibull function was fitted to the data:

$$P(\mathfrak{c}) = 1 - 0.5 \ast \mathrm{e}^{-\left(\frac{\mathfrak{c}}{\mathfrak{a}}\right)^{\beta}} \tag{1}$$

Where P(c) is the probability of correct responses at contrast level c, α is the contrast level that supports threshold performance (82% correct), and β is the slope of the function.

### Cell Mapping and Classification

LGN cells were hand mapped using first a flashlight and then an elongated bar with a sharp contrast profile. Cells were classified as ON-center or OFF-center cells, and Parvocellular (P) or Magnocellular (M) cells, based on their visually driven responses (Norton and Casagrande, 1982; Norton et al., 1988; Xu et al., 2001; Royal et al., 2006; Jiang et al., 2015).

### Neurometric Functions

Basic procedures in computing neurometric functions were similar to those described in previous studies (Barlow et al., 1971; Britten et al., 1992; Purushothaman and Bradley, 2005). For every contrast level, an ROC (Receiver Operating Characteristic) curve was computed (Green and Swets, 1966). Each ROC curve plotted, for all possible signal detection criteria (spikes), the proportion of stimulus-inside-receptive-field trials where the spike count exceeded a certain criterion, against the proportion of stimulus-outside-receptive-field trials that exceeded the same criterion. Next each area-under-ROC curve value was calculated and plotted against its corresponding contrast, a Weibull function ([1], as described above) was fitted, and the neurometric threshold and slope were obtained from the fitted curve. All the LGN cells that could be clearly mapped and maintained long enough to characterize both the psychophysical and neural responses (>150 trials, overall psychophysical performance >65% correct) in the detection task were included in this analysis (overall: n = 89 neurons; monkey 1: n = 61; monkey 2: n = 28). We identified in this dataset 41 ON-center P neurons, 27 OFF-center P neurons, 19 ON-center M neurons, and 2 OFF-center M neurons. We found that the average neurometric threshold (54.4 ± 4.78% contrast, n = 89 neurons, in 0–150 ms integration time windows) was significantly different from the simultaneously measured psychometric threshold (5.76 ± 0.72% contrast, P = 0.000, Wilcoxon signed rank test). The average ratio of neurometric to psychometric threshold was 40.74 ± 10.75, indicating that the average LGN neuron was much less sensitive than the monkey in contrast detection (Jiang et al., 2015).

### Choice Probability

Basic procedures in computing choice probabilities also were similar to those described in previous studies (Britten et al., 1996; Purushothaman and Bradley, 2005). The choice probability for a certain contrast was measured by plotting, as an ROC curve, the proportion of choice-insidereceptive-field trials (i.e., trials in which the monkey saccaded towards the receptive field location) against the proportion of choice-outside-receptive-field trials that exceeded the spike count criteria, and computing the area under that curve. The significance of individual or population choice probabilities was assessed using permutation tests (Britten et al., 1996; Jiang et al., 2015). To accurately estimate choice probability, only neural recordings that met the following criteria were included in this analysis: (1) Behavior ratio (choice-inside/choice-outside-receptive-field trials) > 0.25 and <4; and (2) For every contrast level that was included in the choice probability computation, at least 10 choiceinside and 10 choice-outside-receptive-field trials were recorded. Out of the 89 neurons in the above dataset, 75 (54 P neurons, 21 M neurons) were included in the choice probability analysis according to these criteria. We found in this dataset that, in the absence of any physical stimulus (i.e., 0% contrast, blank trials only), the average choice probability was 0.54 ± 0.01 for LGN P neurons and 0.54 ± 0.01 for LGN M neurons, both above chance according to permutation tests (P neuron: P = 0.015, M neuron: P = 0.033; Jiang et al., 2015).

### Pooling Model

The basic structure of our pooling model was similar to other bottom up pooling models previously proposed to account for the psychophysical threshold and choice probabilities measured during behavioral tasks (Shadlen et al., 1996; Purushothaman and Bradley, 2005; Cohen and Newsome, 2009; Haefner et al., 2013; Liu et al., 2013). Briefly, to simulate a perceptual decision pool of n units, n single neurons were randomly chosen, with replacement, from our entire dataset. To construct a single trial at a given contrast, we simulated each neuron's response by randomly drawing a number from a Gaussian distribution; the mean and variance of this distribution were determined by that neuron's measured response at that contrast level. In each trial, the model made a ''choice'' by comparing the summed activity of the neural pool at the test contrast level to the summed activity of the same neurons at the reference contrast (i.e., blank, 0% contrast trials). This procedure was repeated 50 times (i.e., to simulate 50 trials) for each of the 5 contrast levels, and the simulated ''psychophysical'' performance was recorded as the percentage of correct ''choices'' at each contrast. This performance was fitted with a cumulative Weibull function [1], and threshold and slope parameters were extracted as described above. The choice probability for each simulated neuron was quantified as the covariation between the simulated neural response and the simulated ''psychophysical choice'' at 0% contrast. The parameters for the model included the number of neurons in the pool (n = 1–512 neurons), the integration time window (t = 25–200 ms), the Fano factor (f = 0.25–3.0), the interneuronal noise correlation (r = 0–0.3), and the downstream pooling noise (p = 0–4.0). For each parameter combination, the set of simulations described above (50 trials ∗ 5 contrast levels) was repeated 200 times, each time with a new random sample of n neurons (with replacement), thus giving reliable estimations of the model performance. The overall fitness of the model was evaluated by computing a Goodness-of-Fit (GoF) index:

$$\begin{aligned} \text{GofF} &= (1 - 1/3 \ast |\text{simulated threshold} - \text{measured}| \\ \text{threshold}| / \text{measured threshold} &+ |\text{simulated P choice} \\ \text{probability} &- \text{measured P choice probability} | / \text{measured} \\ \text{P choice probability} &+ |\text{simulated M choice probability} \\ &- \text{measured M choice probability} | / \text{measured M choice} \\ \text{probability} &() \ast 100\% \end{aligned}$$

### Model Parameters

In our simulations we typically used fixed or experimentally measured values (Jiang et al., 2015) for variables such as the Fano factor, the interneural correlation, and the pooling noise (with the exception of **Figure 1**, where these parameters were systematically varied to probe the basic properties of the pooling model). We consider our choices for these parameter values to be neurobiologically realistic and meaningful for the following reasons: (1) The Fano factor: Recordings in anesthetized as well as alert animals have reported significant variabilities in the responses of single cortical neurons, with the Fano factor (response variance/mean) averaging 1.0–3.0 (Tolhurst et al., 1983; McAdams and Maunsell, 1999; Oram et al., 1999; Gu et al., 2007). The Fano factor of subcortical visual neurons, however, is relatively low (i.e., <1.0). This is true for retinal ganglion cells (Levine et al., 1992; Berry et al., 1997; Reich et al., 1997) as well as LGN cells (Kara et al., 2000). The Fano factors measured in our detection task (Jiang et al., 2015) and used in our models (0.8–1.4, depending on integration time) were in agreement these previously reported measurements (2) Interneuronal correlation: In sensory cortex, interneuronal correlations between pairs of nearby neurons are typically weak but positive (∼0.1–0.2; Averbeck et al., 2006). For the LGN P-P and M-M neuron pairs, because convergent feedforward, divergent feedforward, and lateral connections are sparser than those in the cortex (Casagrande and Xu, 2004; Nassi and Callaway, 2009), it is not surprising that we found an average interneuronal correlation (0.028) that was not significantly different from 0.0. The interneuronal correlation for a P-M neuron pair is very likely to be even smaller, as the P and M pathways receive different retinal inputs, remain segregated in different layers of the LGN (Casagrande and Norton, 1991; Nassi and Callaway, 2009), and retain separate feedback loops with V1 (Ichida and Casagrande, 2002; Briggs and Usrey, 2009, 2011; Ichida et al., 2014). In our model the P-M correlation was fixed at 0.01, but our simulations could always approach >99% GoF at some parameter combinations, given any P-M correlation values between 0.0–0.05 (data not shown). (3) Pooling noise: The downstream pooling noise can be thought of as the average Fano factor of the cortical neurons onto which LGN neurons converge (Shadlen et al., 1996). In our simulations this number was fixed at 2.0, which is the average estimation of the Fano factor in cortex (see above). The success of our simulations (i.e., approaching >99% GoF at some parameter combinations), however, did not . depend on this assumption. Similar model results could be obtained by simply assuming a true Poisson distribution for all downstream neurons (i.e., pooling noise = 1.0).

### Uniform Pooling and Alternative Pooling Strategies

In this paper, we examined uniform pooling models as well as several alternative weighted pooling models. All of these models shared the same overall structure:

$$X\_{\text{pooled}} = \sum\_{i=1}^{n} w\_i x\_i \tag{3}$$

Where Xpooled is the summed activity of the perceptual decision pool, xi is the response of a single neuron, and wi is the readout weight assigned to this neuron. Within this structure, the uniform pooling model simply assigned equal weights (i.e., wi = 1.0) to all neurons in the decision pool, regardless of their sensitivities (Shadlen et al., 1996). An alternative weighted pooling strategy, in contrast, calculated the sensitivity of each individual neuron and assigned weights accordingly. Depending on how this neural sensitivity was quantified, there were three main categories of weighted pooling schemes: (1) Amplitudeper-trial (amp/trial) weighted scheme, where every neuron was weighted according to its response amplitude in every trial (wi ∝ xi), with the neuron with the highest spike rate carrying a weight of 1.0; (2) Mean amplitude (mean amp) weighted scheme, where every neuron was weighted according to its average response amplitude (wi ∝ xi) at high contrast (80–99%), with the neuron with the highest average spike rate carrying a weight of 1.0; and (3) D-prime weighted scheme, where every neuron was weighted by its d-prime value (wi ∝ d 0 i ) at high contrast (80–99%), with the neuron with the greatest d-prime carrying a weight of 1.0. Here d-prime was defined as:

$$d\_i^{'} = (\overline{x}\_i - \overline{x}\overline{0}\_i) \bigg/ \sqrt{\frac{s\_i^2 + s\mathbb{O}\_i^2}{2}} \tag{4}$$

Where xi is the neuron's mean response amplitude at high contrast (80–99%), x0i is its mean response amplitude at reference contrast (i.e., blank, 0% contrast), and si and s0i represent the corresponding standard deviations. Because of their potential deviations from normality, for both weight and d-prime distributions we reported medians as well as means. Additionally, to characterize the spread of a distribution, we reported the interquartile range:

$$\text{Interquartile Range } (IQR) \; = \; Q\_3 - Q\_1 \tag{5}$$

Where Q1 is the 1st quartile (i.e., 25% percentile), and Q3 is the 3rd quartile (i.e., 75% percentile) in the range. To characterize the skewness of a distribution, we reported the skewness index:

$$\text{Skewness index (SI)} = \frac{\frac{1}{n} \sum\_{i=1}^{n} \left(\chi\_i - \overline{\chi}\right)^3}{\left(\sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left(\chi\_i - \overline{\chi}\right)^2}\right)^3} \tag{6}$$

Where n is the sample size, and x is the mean of the sample distribution.

noise was assumed to be 0; integration time was 0–150 ms. As in (B), each simulation line moved from the upper right corner to the lower left corner as more neurons were added in the pool. Legends as in (C,D). (G,H) Increasing pooling noise increased simulated threshold and decreased choice probability, for both P neurons (G) and M neurons (H). The Fano factor was assumed to be 1.03 (measured value); interneuronal correlation was assumed to be 0.028 (measured value); integration time was 0–150 ms. As in (B), each simulation line moved from the upper right corner to the lower left corner as more neurons were added in the pool. Legends as in (C,D). (A,B) and (E–H) were adapted from Jiang et al. (2015).

### Results

As previously reported (Jiang et al., 2015), we found that single LGN P and M neurons, although not as sensitive as the monkeys in detecting contrast, were significantly correlated with the behavioral choices of the monkeys during a 2AFC contrast detection task. Based on these experimental data, we report in this paper a series of modeling results in the following order: first we describe the basic parameters of the uniform pooling model and its performance in different time frames, then we compare several alternative weighted pooling schemes to the uniform pooling model, and finally we zoom in on one of the best performing weighted schemes and examine its structure in detail.

### Uniform Pooling Model: Parameters

The uniform pooling model we built was similar to a number of previous models used to account for psychophysical performance and choice probability measurements based on sensory neural responses (Shadlen et al., 1996; Purushothaman and Bradley, 2005; Cohen and Newsome, 2009; Haefner et al., 2013; Liu et al., 2013). Inputs to the model were single LGN P and M neural responses at different contrasts. Outputs were the simulated ''psychophysical'' threshold and the simulated choice probabilities for individual neurons. In agreement with Shadlen et al. (1996), this uniform pooling model behaved predictably when certain model parameters were changed. Specifically, increasing the number of neurons in the pool (n) decreased the psychophysical threshold and choice probability values (**Figures 1A,B**). Increasing the Fano factor increased the simulated threshold but maintained the same choice probability values (**Figures 1C,D**). Increasing the interneuronal correlation increased the threshold as well as choice probability values (**Figures 1E,F**). Increasing the downstream pooling noise increased the threshold and decreased choice probability values (**Figures 1G,H**). For the simulation results reported in the following sections, the Fano factor and interneuronal noise correlation were fixed at experimentally measured values and the pooling noise, which could be considered as the Fano factor of downstream neurons onto which LGN neurons converge, was assumed to be 2.0 (see ''Material and Methods'' Section).

### Uniform Pooling Model: Performance

The first pooling scheme we investigated was the uniform pooling model, in which the responses of all neurons were weighted equally and summed up to form perceptual decisions. In this model, the simulated psychometric threshold consistently decreased as: (1) more neurons (n) were added into the pool; and (2) the integration time window (t) was extended. Specifically, we found that: (1) at extremely short intervals (25 ms), incorporating a large number of neurons from both the P and M populations (n = 512 P neurons, 512 M neurons) still failed to achieve great psychophysical sensitivities (i.e., threshold <10% contrast), but incorporating a large number of M neurons (n = 256–512) rather than P neurons was more beneficial to model performance (**Figure 2A**); (2) at relatively brief intervals (50 ms), preferably incorporating a large number of M neurons (n = 256–512) rather than P neurons produced good psychophysical performance (i.e., threshold <10% contrast; **Figure 2B**); (3) at medium to long intervals (75–200 ms), a wider range of M/P neuron combinations (n = 64–512) yielded good model performance (**Figures 2C–F**).

Next we analyzed the simulated choice probabilities for the P and M populations and compared them to the measured choice probability distributions (Jiang et al., 2015). Briefly, in the 2AFC contrast detection task, we found that in the absence of any physical stimulus (i.e., 0% contrast, blank trials only), the average choice probability was 0.54 ± 0.01 for LGN P neurons and 0.54 ± 0.01 for LGN M neurons, both above chance according to permutation tests (P neuron: P = 0.015, M neuron: P = 0.033). In the uniform pooling model, the simulated choice probability distributions for P and M neurons (n = 512 neurons, t = 0–150 ms, P choice probability = 0.54 ± 0.00, M choice probability = 0.54 ± 0.00) resembled their experimentally measured counterparts (P > 0.05, permutation tests; **Figures 3A,B**). These choice probability patterns remained unchanged throughout the 200 ms stimulus presentation time (n = 512 neurons, P > 0.05, 1-way ANOVAs; **Figures 3C,D**).

To evaluate the overall performance of this model, a Goodness-of-Fit (GoF) index (see equation [2]) was reported for each (n, t) parameter combination. A GoF (ranging from 0–100%) reflected three factors equally: (1) how close the simulated ''psychophysical'' threshold approached the measured psychophysical threshold; (2) how close the simulated P population choice probability approached the measured average choice probability for P neurons; and (3) how close the simulated M population choice probability approached the measured average choice probability for M neurons. A GoF of 100% indicated that our simulation perfectly matched the observed psychometric threshold and the choice probabilities for both types of neurons. By changing the duration of the integration window (25–200 ms), we found that: (1) at extremely short intervals (25 ms), even incorporating a large number of neurons from both groups (n = 512 P neurons, 512 M neurons) still failed to reproduce the observed threshold and choice probabilities (**Figure 4A**); (2) in 50 ms, preferably incorporating a large number of M neurons (n = 256–512) could explain the observed threshold and choice probabilities (**Figure 4B**); (3) in 75 ms, incorporating a large number of either P or M neurons (n = 128–512) could achieve good overall model performance (**Figure 4C**); and (4) at medium to long intervals (100–200 ms), a smaller number of P and M neurons were needed (n = 32–128) to achieve good model performance, but further increasing the number of neurons resulted in a decrease in model performance (**Figures 4D–F**).

### Alternative Pooling Schemes

In this section we examine several alternative pooling schemes where, instead of assigning the same weight to all neurons, each individual neuron was weighted differentially based on its response rate or sensitivity. We investigated three main categories of weighted pooling schemes, namely the

amplitude-per-trial (amp/trial) weighted, the mean amplitude (mean amp) weighted, and the d-prime weighted schemes (see ''Material and Methods'' Section). First we compared the simulated psychometric thresholds and found that for all pooling schemes the average psychometric threshold decreased with time (n = 1–512 neurons, F = 920.22, P = 0.00, 2 way ANOVA main effect for time). Furthermore, there was a significant difference in psychometric thresholds among different pooling schemes (F = 54.82, P = 0.00, 2-way ANOVA main effect for pooling strategy), and this difference changed across time (F = 8.43, P = 0.00, 2-way ANOVA interaction effect; **Figure 5A**). Next we examined whether these alternative pooling strategies improved the sensitivity of the model when compared to the uniform pooling strategy. Here the mean amplitude weighted and the d-prime weighted models could be further divided into two subcategories, respectively, depending on whether P and M neurons were weighted separately or together in the model. Among all of these alternative pooling schemes, we found that only the d-prime weighted schemes consistently improved the psychophysical performance when compared with the uniform pooling scheme (mean difference, d-prime 1 = −4.88 ± 0.42% contrast;

probability distributions for P (magenta) and M (green) neurons in 200 simulations (n = 512 neurons, t = 0–150 ms). (B) Choice probability distributions for P (magenta) and M (green) neurons in 200 simulations (n = 512 neurons, t = 0–150 ms). Arrow: mean choice probability; solid line: choice probability = 0.5. (C) Cumulative choice probability distributions for P (magenta) neurons in different integration time windows (n = 512 neurons, 200 simulations for each integration time). (D) Cumulative choice probability distributions for M (green) neurons in different integration time windows (n = 512 neurons, 200 simulations for each integration time).

mean difference, d-prime 2 = −4.67 ± 0.43% contrast; P < 0.05, Tukey's HSD tests for multiple comparisons). The mean amplitude weighted schemes and the amplitude per trial weighted scheme all failed to perform as well as the uniform pooling scheme in terms of the threshold (mean difference, mean amp 1 = 2.59 ± 0.21% contrast; mean difference, mean amp 2 = 3.31 ± 0.23% contrast; mean difference, amp/trial = 7.51 ± 0.28% contrast; P < 0.05, Tukey's HSD tests for multiple comparisons). Additionally, the two subtypes of mean amplitude weighted models did not differ from each other in terms of their simulated thresholds, and the two subtypes of d-prime weighted models did not differ from each other either (P > 0.05, Tukey's HSD tests for multiple comparisons; **Figure 5B**). Finally, the minimal psychophysical threshold achieved by the model also decreased with time in all pooling schemes and plateaued at around 50–75 ms after stimulus onset (n = 1–512 neurons, minimal threshold = 2–3% contrast; **Figure 5C**). Taken together, **Figures 5A–C** demonstrated that the d-prime weighted pooling strategies were the most optimal in terms of the simulated psychophysical performance, and this advantage over other pooling strategies was the most apparent in short integration time windows (25–50 ms).

Next, we compared the simulated choice probabilities in these different pooling schemes. First, in a fixed time window (t = 0–150 ms), the overall choice probability distributions for simulated P or M neurons did not differ significantly among the pooling schemes (n = 512 neurons, P > 0.05, 1-way ANOVAs; **Figures 5D,E**). Second, these choice probability distributions did not change with time in the case of the amplitude per trial and mean amplitude weighted models (n = 512 neurons, P > 0.05, 2-way ANOVAs main effect for time; **Figures 5F,G**). In the d-prime weighted pooling model, however, the choice probability distributions did shift significantly as the integration time window was extended (n = 512 neurons, F = 4.86, P = 0.00, 2-way ANOVA main effect for time; **Figure 5H**), corresponding well to the temporal dynamics of choice probability that were experimentally measured in LGN P and M neurons (see Figure 7 of Jiang et al., 2015).

Finally, we investigated the overall performance of different pooling strategies by comparing their GoF indices in different time windows. We found that for all pooling schemes the overall