# FROM CELL PHYSIOLOGY TO EMERGING BRAIN FUNCTIONS

EDITED BY : Philippe Isope, Lisa Mapelli, Marylka Yoe Uusisaari and Alexander Tang PUBLISHED IN : Frontiers in Cellular Neuroscience

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88966-000-1 DOI 10.3389/978-2-88966-000-1

#### 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>
# FROM CELL PHYSIOLOGY TO EMERGING BRAIN FUNCTIONS

Topic Editors:

Philippe Isope, Centre National de la Recherche Scientifique (CNRS), France Lisa Mapelli, University of Pavia, Italy Marylka Yoe Uusisaari, Okinawa Institute of Science and Technology Graduate University, Japan Alexander Tang, University of Western Australia, Australia

Citation: Isope, P., Mapelli, L., Uusisaari, M. Y., Tang, A., eds. (2020). From Cell Physiology to Emerging Brain Functions. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88966-000-1

# Table of Contents

Elena Laura Georgescu, Ioana Antoaneta Georgescu, Carmen Denise Mihaela Zahiu, Alexandru Răzvan Şteopoaie, Vlad Petru Morozan, Adrian Ştefan Pană, Ana-Maria Zăgrean and Daniela Popa

*112 How Do Cells of the Oligodendrocyte Lineage Affect Neuronal Circuits to Influence Motor Function, Memory and Mood?*

Renee E. Pepper, Kimberley A. Pitman, Carlie L. Cullen and Kaylene M. Young

*126 PDI Knockdown Inhibits Seizure Activity in Acute Seizure and Chronic Epilepsy Rat Models via S-Nitrosylation-Independent Thiolation on NMDA Receptor*

A. Ran Jeon and Ji-Eun Kim

Daniel Aharoni and Tycho M. Hoogland

*195 Spontaneous and Acetylcholine Evoked Calcium Transients in the Developing Mouse Utricle*

Holly A. Holman, Lauren A. Poppi, Micah Frerck and Richard D. Rabbitt

*213 Enriched Environment Shortens the Duration of Action Potentials in Cerebellar Granule Cells*

Abdelmoneim Eshra, Petra Hirrlinger and Stefan Hallermann

# Maturation of GABAergic Transmission in Cerebellar Purkinje Cells Is Sex Dependent and Altered in the Valproate Model of Autism

Sébastien Roux 1 , Ann Lohof 2 , Yehezkel Ben-Ari 3 , Bernard Poulain1 and Jean-Louis Bossu1 \*

1 Institut des Neurosciences Cellulaires et Intégratives (INCI)—CNRS, UPR 3212, Strasbourg, France, 2Sorbonne Université, CNRS UMR 8256, Biological Adaptation and Ageing, Paris, France, 3Neurochlore, Ben-Ari Institute of Neuroarcheology, Campus Scientifique de Luminy, Aix Marseille Université, Marseille, France

Brain development is accompanied by a shift in gamma-aminobutyric acid (GABA) response from depolarizing-excitatory to hyperpolarizing-inhibitory, due to a reduction of intracellular chloride concentration. This sequence is delayed in Autism Spectrum Disorders (ASD). We now report a similar alteration of this shift in the cerebellum, a structure implicated in ASD. Using single GABAA receptor channel recordings in cerebellar Purkinje cells (PCs), we found two conductance levels (18 and 10 pS), the former being dominant in newborns and the latter in young-adults. This conductance shift and the depolarizing/excitatory to hyperpolarizing/inhibitory GABA shift occurred 4 days later in females than males. Our data support a sex-dependent developmental shift of GABA conductance and chloride gradient, leading to different developmental timing in males and females. Because these developmental sequences are altered in ASD, this study further stresses the importance of developmental timing in pathological neurodevelopment.

Keywords: cerebellum, Purkinje neuron, chloride gradient, autism, sexual dimorphism, GABA**A** channels

# INTRODUCTION

Gamma-aminobutyric acid (GABA) is the main inhibitory neurotransmitter in the central nervous system. By opening GABAA receptor channels and via chloride flow down its electrochemical gradient, GABA hyperpolarizes the membrane potential and thus reduces neuronal excitability in the mature nervous system. Intracellular chloride concentration (Cl−)I is regulated by cationchloride co-transporters that determine the strength of GABAergic inhibition (Rivera et al., 1999; Watanabe and Fukuda, 2015; Raimondo et al., 2017). If the intracellular chloride concentration is high enough that the equilibrium potential for chloride is positive compared to the resting membrane potential, GABAA receptor activation can depolarize the cell above the action potential threshold, acting like an excitatory neurotransmitter (Ben-Ari et al., 1989; Ben-Ari, 2014). This regulation of (Cl−)I has been reported in a wide range of animal species and brain structures; it is affected by many factors including sex steroids, BDNF and IGF-1 signaling (Galanopoulou, 2008; Tsutsui et al., 2011; Waddell and McCarthy, 2012; Ben-Ari, 2014; Watanabe and Fukuda, 2015). Immature neurons have lower (Cl−)I than neurons of young adults, leading to depolarizing GABA actions which both activate voltage-gated calcium currents and allow NMDA receptor activation, underlying the well-known trophic actions of GABA in development (Ben-Ari, 2002, 2014; Witte et al., 2014; Raimondo et al., 2017). GABAergic signals also generate action potentials

#### Edited by:

Enrico Cherubini, Scuola Internazionale Superiore di Studi Avanzati (SISSA), Italy

#### Reviewed by:

Tommaso Pizzorusso, Consiglio Nazionale delle Ricerche (CNR), Italy Andrea Barberis, Fondazione Istituto Italiano di Technologia, Italy

\*Correspondence:

Jean-Louis Bossu<EMAIL_ADDRESS>
Received: 04 May 2018 Accepted: 16 July 2018 Published: 30 July 2018

#### Citation:

Roux S, Lohof A, Ben-Ari Y, Poulain B and Bossu J-L (2018) Maturation of GABAergic Transmission in Cerebellar Purkinje Cells Is Sex Dependent and Altered in the Valproate Model of Autism. Front. Cell. Neurosci. 12:232. doi: 10.3389/fncel.2018.00232 in many immature neurons, and the switch between excitatory and inhibitory GABA effects occurs at specific developmental time-points, depending on the neural structure and animal species. For example, the GABA ''switch'' occurs between E15.5 and E 17.5 in mouse spinal motoneurons (Allain et al., 2011), around P5 in the ventral horn of the rat spinal cord (Stein et al., 2004) and in ganglion and amacrine cells of the retina (Zhang et al., 2006), and from P8–P10 in mouse hippocampal pyramidal cells (Ben-Ari et al., 1989; Tyzio et al., 2007). Changes in the timing of the GABA switch, due to genetic mutations and/or environmental insults, are thought to contribute to neurodevelopmental disorders including autism (Ben-Ari et al., 1989; Ben-Ari, 2017).

GABAA receptors are composed of 2 α subunits, 2 β subunits and 1Υ subunit, which together form a central ion pore. Excluding splice variants and point mutations, 6 α isoforms, 3 β isoforms and 3 Υ isoforms have been characterized, as well as some minor subunits (δ, ρ, ε, θ, π and ρ). In theory, therefore, a very large number of GABAA receptor types may be found even in a single cell. The major adult isoform is generally accepted to be composed of α1, β2 and Υ2 subunits. Whereas some subunits have a broad expression throughout the central nervous system, other subunits show a restricted cellular and subcellular localization (Sigel and Steinmann, 2012). The expression of GABAA receptor subunits is developmentally regulated, notably with a developmental switch in α subunit expression associated with a slow-to-fast shift in the kinetics of GABA-mediated inhibitory postsynaptic potentials (Laurie et al., 1992; Fritschy et al., 1994). Thus, in thalamic reticular neurons, a postnatal switch in GABAA receptor subunits from α5 to α3 is believed to play a role in the early development of the circuit (Pangratz-Fuehrer et al., 2016). Experiments using subunit expression in HEK cells or mouse fibroblasts show that single GABA channels currents have distinct opening conductance levels depending upon subunit composition (Mortensen and Smart, 2006).

Studies on developmental GABA shifts during normal development and in relation to autism spectrum disorders (ASD) have been done primarily in cortical structures. Yet abnormalities of the cerebellum and of cerebellar Purkinje cells (PCs) are frequently described in post mortem studies of humans with autism (Fatemi et al., 2012) as well as in the rodent valproate model of autism (Ingram et al., 2000). In addition, PCs undergo considerable post-natal development, including regression of climbing fiber multi-innervation and dramatic dendritic arbor expansion (Dusart and Flamant, 2012), allowing investigation of these important developmental processes after birth.

We have now measured the changes in somatic GABAA channel properties and chloride gradient in PCs, comparing these changes in male and female mice, since sex differences in the GABA shift have been reported (Galanopoulou, 2008). We report that in normal and valproate-model male and female mice, the dominant GABAA receptor channel shifts from high-conductance in the newborn to low-conductance in juvenile and adult mice. This shift in conductance parallels the GABA switch from depolarizing to hyperpolarizing, as it is sex dependent, being delayed in naïve females compared to males. Furthermore, the GABA switch is delayed in mice of both sexes after prenatal exposure to valproate, supporting a role of cerebellar dysfunction in the pathology of autism.

# MATERIALS AND METHODS

## Mice

Pregnant C57 mice were injected intraperitoneally with 600 mg valproate sodium salt (n = 18) or saline (n = 15) at embryonic day 12.5 (Roullet et al., 2013). Either control pups or pups from valproate-treated dams were used to prepare acute cerebellar slices for patch-clamp experiments at ages between postnatal days 5 (P5) and 45 (P45). For the study of climbing fiber synapse elimination, either C57 or Swiss pups (because of large number of pups per litter) were used; no differences were found in synapse elimination between the two strains of mice. All procedures followed guidelines established by le Comité National d'ethique pour les Sciences de la Vie et de la Santé (EU Council Directive 2010/63/EU) and were approved by Institutional Animal Care and Use Committees (CREMAS, Comité Régional d'Ethique en experimentation animale de Strasbourg).

#### Slice Preparation

Standard procedures were used to prepare 250-µm or 300-µm parasagittal slices from control or valproate-treated mice at P5–P45 following a protocol approved by the European and French guidelines on animal experimentation established by le Comité National d'ethique pour les Sciences de la Vie et de la Santé (EU Council Directive 2010/63/EU) and were approved by Institutional Animal Care and Use Committees (CREMAS, Comité Régional d'Ethique en experimentation animale de Strasbourg). Briefly, mice were killed by decapitation under isoflurane anesthesia. Brains were dissected in ice-cold artificial cerebrospinal fluid (ACSF) and sliced with a vibratome (Leica VT1200S) at 4◦C. Slices were maintained for 30 min at 32◦C in an interface chamber containing ACSF equilibrated with 95% O2, 5% CO2 and containing (in mM): NaCl 124, KCl 2.7, CaCl2 2, MgCl2 1.3, NaHCO3 26, NaH2PO4 0.4, glucose 10, ascorbate 4, then for at least 1 h at room temperature before being transferred to a superfusing recording chamber.

# Electrophysiological Recordings

Slices were transferred to a recording chamber on an upright microscope. The recording chamber was continuously perfused at room temperature with bath solution containing: (mM) NaCl 124, KCl 2.7, CaCl2 2, MgCl2 1.3, NaHCO3 26, NaH2PO4 0.4, glucose 10, pH 7.4, equilibrated with 95% O2, 5% CO2. For cell-attached recordings (in order to stabilized the resting membrane potential) the bath solution contained tetrodotoxin (TTX) 10−5 M and NBQX 10−5 M. In some experiments isoguvacine or NBQX were applied to the bathing fluid at a concentration of 10−5 M. For experiments recording climbing fiber currents, the ACSF contained 10−4 M picrotoxin.

Most electrophysiological experiments were performed on visually identified PCs using the patch-clamp technique in the cell-attached configuration. Electrodes were filled for single channels recordings with the following solution (mM): KCl 110, NaCl 2, MgCl2 2, CaCl2 2, HEPES 10, tetra-ethyl-ammoniumchloride (TEA) 20, TTX 10−3 , CsCl2, 4-aminopyridine (4AP) 1, BaCl2 1, Isoguvacine or Muscimol 10−5 , pH 7.4 ; and for spiking activity with the following solution (mM): NaCl 124, KCl 2.7, CaCl2 2, MgCl2 1.3, NaHCO3 26, NaH2PO4 0.4, glucose 10, equilibrated at pH 7.4 with 95%O2, 5% CO2.

For experiments recording climbing fiber currents, patch pipettes were filled with a solution containing (mM): Cs-D-gluconate 120, biocytin 13, 10 HEPES, BAPTA 10, TEACl 3, Na2ATP 2, MgATP 2, NaGTP 0.2, pH 7.3, 290–300 mOsm. Climbing fiber currents were elicited by stimulation in the internal granular layer with a saline-filled glass pipette.

Signals were recorded and filtered at 5 kHz using an Axopatch 200A amplifier (Axon Instrument). Current and voltage signal were digitized at 50 kHz using a Digidata 1322A (Axon Instruments) prior to being recorded directly using Clampex (10.2) software. Analysis was performed off-line using Clampfit (10.2) software. Data were filtered before analysis with a cut off frequency of 1.5 or 1 kHz.

#### Graphs, Fitting Procedures and Statistics

Sigma plot 12.5 software was used for graphic representations of the data, fitting procedures and statistical analysis. For channels conductance analysis the normality Shapiro-Wilk test, and equal variance test were used before running an un-paired Student's t-test. Data were considered statistically significant when P < 0.05. ∗ is used for P values between 0.05 and 0.01; ∗∗ for P values between P > 0.01 and 0.001. ∗∗∗ is used for P values < 0.001. The quality of the fit was determined using the prediction error.

# RESULTS

## Two GABAA Conductances in Neonatal PCs

The single-channel properties of GABAA receptors have been extensively studied in neuronal cultures and brain slices. In cerebellar granule cells, three main conductance (28, 17, 12 pS) have been characterized and are attributed to distinct GABAA receptor subtypes (Brickley et al., 1999). We analyzed currents activated by isoguvacine or muscimol, two GABAA receptor agonists, recorded from 108 cell-attached patches from PCs in acute cerebellar slices from male and female mice, between P5 and P45. In 10 patches, two conductance levels were recorded simultaneously (**Figure 1**). **Figure 1A1** illustrates current traces showing two channels subtypes with a slope conductance of 19 pS and 7.3 pS, reversing at the same potential, indicating a similar ionic selectivity (**Figure 1A3**). The amplitude distribution in **Figure 1A2** shows that the dominant channel conductance in this patch is level 1, and that it is not the result of simultaneous openings of two level 2 channels. The mean conductance slope value is 17.7 ± 0.9 pS for level 1 and 8.1 ± 0.5 pS for level 2 (**Figure 1B**, n = 10); these are significantly different (P < 0.001).

The histogram distribution of the slope conductance of level 1 (in red) and 2 (in blue) channels are illustrated in **Figure 1C**. For both channel types the amplitude distribution is normal and fitted with a Gaussian function but with distinct peak values, 17 pS for level 1 and 7.5 pS for level 2.

# GABA Conductance Switch During PC Development

A switch in the dominant GABAA receptor channel conductance occurs during development, from primarily level 1 in immature PCs (P5–12) to primarily level 2 in mature (P26–45) PCs in both sexes. Two recordings from PCs (female mice) illustrating this switch are shown in **Figure 2A1**. At P6, the slope conductance of the dominant GABAA receptor channel is 18 pS, whereas at P24 it is 10.8 pS (**Figure 2A2**).

**Figure 2B** summarizes and compares the conductances of the dominant channels recorded in males (left panel) at P5–12 (newborn) and at P28–45 (young adult) and in females (right panel) at P6–12 and P26–36. In males the mean conductance of the dominant GABAA channel switch from 17.3 ± 1.6 pS (n = 14) at P5–12 to 10.5 pS ± 1.7 pS (n = 15, p < 0.001) at P28–45; and in females from 18.5 ± 3.0 pS (n = 11) at P6–12 to 9.8 ± 8 pS (n = 6, p < 0.001) at P26–36. In newborn mice (both sexes), the mean conductance of the dominant GABAA channel is not significantly different from the level 1 conductance illustrated in **Figure 1**; similarly, the mean conductance of the dominant GABAA channel in young-adult mice (both sexes) is not significantly different from the level 2 conductance illustrated in **Figure 1**. Thus GABAA channels in PCs switch from level 1 to level 2 during post-natal development.

#### Sexual Dimorphism of the Chloride Gradient Shift

We then determined the reversal potential of the dominant GABAA receptor channel relative to the resting membrane potential (RP) during PC development in males and females. **Figure 3** illustrates current recordings obtained in PCs from a P16 male (**Figure 3A**) and a P16 female (**Figure 3B**). The reversal potential of the GABAA receptor channel is determined by the linear regression used to fit the current/voltage (I/V) relationship of a single GABAA receptor channel. At P16, the reversal potential is negative to RP (−14 mV) in males (**Figure 3A**, bottom panel), but positive to RP (+4 mV) in females (**Figure 3B**, bottom panel).

The changes in GABA reversal potential relative to RP in cerebellar PCs from males and females is illustrated in **Figure 3C**. The curve in males can be fitted with a sigmoidal function (four parameters, X0 14.4 days). Between P5–P11 the reversal potential is stable, around +10 mV relative to RP; the reversal potential then rapidly reaches a negative value around P14 and stabilizes at P26 around −10 mV relative to RP (−10.3 ± 5.0 mV, n = 10 at P ≥ 26, see bottom graph in inset). Therefore, GABA exerts depolarizing effects in male PCs prior to P15 and hyperpolarizing effects thereafter.

FIGURE 1 | Two levels of type A gamma-aminobutyric acid (GABAA) receptor channel conductances can be recorded from the same Purkinje cell (PC) membrane patch. (A1) Current traces recorded in cell-attached conditions at a potential of 80 mV relative to resting potential (RP). Two amplitude levels of GABAA receptor channels are seen (level 1 in red, level 2 in blue). (A2) Data from the same patch as in panel (A1) is illustrated in a current amplitude histogram. The distribution is fitted by two Gaussians, with peaks at 1.65 pA and 0.74 pA. (A3) The relation between single channel current amplitude and the potential relative to RP for the two conductance levels. Linear regressions give slope conductances of 19 pS (red) and 7.3 pS (blue). The reversal potential for the two conducting levels is −11 mV relative to RP. (B) Bar graphs of level 1 (red) and level 2 (blue) conductances recorded in 7 PC membrane patches. The mean value of level 1 is 17.7 ± 0.9 pA and differs significantly (∗∗∗p < 0.001) from the mean value of level 2 (8.1 ± 0.5 pA). (C) Histogram distribution of the slope conductance of level 1 (red) and level 2 (blue). Each conductance distribution follows a normal distribution fitted by a Gaussian with a peak value of 17 pS for level 1 and 7.5 pS for level 2.

With the exception of the earliest post-natal days, the developmental curve of the reversal potential in female PCs also follows a sigmoid pattern from P8 to P21 but with a systematically more depolarized value compared to males (**Figure 3C**, red symbols). Although the reversal potential is not significantly more negative between P5 and P7 in females than in males (+3.5 mV vs. +8.7 mV), it then becomes significantly (p < 0.001) more positive between P8 and P12 (+15.8 ± 3.6 mV, n = 11) than in males (+9.1 ± 4.1 mV, n = 8, see top graph in inset). From this higher plateau value, the reversal potential change over time is best fitted by a sigmoid curve (four parameters, X0 18.7 days) and decreases to −8.5 ± 3.5 mV (n = 8), at P26 and older (see bottom graph in inset). These observations show that GABA exerts more depolarizing actions in female PCs than in male PCs during a period of nearly 2 weeks (P8 to P21). In addition, the GABA switch is delayed in females by 4 days (P18 vs. P14 in males), such that from P15–19 the reversal potential is negative relative to RP in males (−15.9 ± 7.4 mV, n = 17) but still positive in females (+7.6 ± 5 mV, n = 10; see top graph in inset).

An important step in cerebellar synaptic maturation during the early postnatal period is the refinement of climbing fiber synapses on PCs, from a multiple climbing fiber innervation to mono-innervation (Hashimoto and Kano, 2003). To determine possible links between GABA signaling and circuit refinement in terms of climbing fiber synapse elimination, we measured

the number of climbing fiber synapses on each PC from males and females at different ages (**Figure 3D**). The progress of synapse elimination in both sexes was very similar, indicating that sex-dependent differences in GABAergic signaling does not disrupt this fundamental developmental process.

# The Developmental Switch of the Chloride Gradient Is Delayed in the Valproate Mouse Model of Autism

In the valproate rat model of autism, the developmental regulation of the chloride gradient, which determines the effects of GABAA receptor activation on neuronal excitability, is disrupted in hippocampal neurons. As a result, the developmental excitatory-inhibitory switch of GABAergic effects is abolished in hippocampal neurons from valproate rats, with depolarizing actions from birth onwards (Tyzio et al., 2014). We used this same model in mice to analyze the developmental changes in GABAA receptor channel conductance and the chloride gradient. We recorded 57 PC-attached recordings of GABAA receptor channels from 35 males and 22 females from valproate-treated dams, at different postnatal ages. As in control mice, recordings from valproate mice revealed two conductance levels (level 1 mean value 15.8 ± 0.4 pS and level 2 mean value 7.8 ± 0.5 pS, p < 0.001) that could occasionally be recorded from the same patch (five patches, **Figure 4A**). Furthermore,

FIGURE 3 | Developmental evolution of the GABAA receptor channel reversal potential in males and females and climbing fiber synapse elimination. (A) Recordings from a P16 male. Current traces show the GABAA receptor channel recorded from a PC membrane patch at three potentials as indicated on top of each trace. The I/V curve below, a linear regression of the data from this channel, gives a conductance of 15.2 pS and a reversal potential with an extrapolated value of −14 mV relative to RP. (B) Recordings from a P16 female. Current traces show the GABAA receptor channel recorded from a PC membrane patch at three potentials as indicated on top of each trace. The I/V curve below, a linear regression of the data from this channel, gives a conductance of 17.5 pS and a reversal potential with an extrapolated value of + 4 mV relative to RP. Note that at the same post-natal age the GABAA receptor channel reverses at a negative potential relative to RP in the male, but in the female reversal occurs at a positive potential relative to RP. (C) Left graph, evolution of the mean reversal potential of the dominant GABAA receptor channel as a function of the age in males (black) and females (red). The data were best fit (R 0.98) by a sigmoidal hill function with four parameters (f = Y0 + a/(1 + exp(−(X − X0)/b)) between P5 and P48 in males (black line, a (max) = 22, b (slope) = −0.31, X0 = 14.4 and Y0 (min) = −13.5) and in between P8 and P36 in females (red line, a = 24.2, b = −1.25, X0 = 18.5, Y0 = −8.3). The horizontal line at 0 mV shows the depolarizing/hyperpolarizing switch. The bar graphs in inset compare the mean values of the reversal potential relative to RP in males (black) and in females (red) at P8–12, P15–19 and P26–36. A significant difference between males and females (∗∗∗p < 0.001) was found at P8–12 and P15–19. (D) Right graph: regression of climbing fiber multi-innervation of Purkinje cells (PCs) is similar in males and females, indicating that this basic circuit maturation during this developmental period is not sex-dependent.

as in control mice, PCs from valproate-treated male or female newborn mice show a dominant level 1 conductance channel, and PCs from juvenile and young adult PCs showed a dominant level 2 conductance channel. In males, the mean conductance switches from 17.1 ± 1.6 pS (n = 13) at P6–10, to 13.1 pS ± 1.8 pS

(n = 5) at P24–30 (**Figure 4B**, p = 0.002). In females, the mean conductance switches from 17.7 ± 3.2 pS (n = 5) at P7–10, to 10.9 pS ± 1.6 pS (n = 7) at P23–30 (**Figure 4C**, p < 0.001). These conductance levels are not different between the sexes nor are they different from conductance levels 1 and 2 (**Figures 1B**, **4A**)

X0 = 23 and Y0 = −5). The red curve is a sigmoidal fit of the evolution of the reversal potential in control females and the red open dots the mean value in control females not included in the fit by the red line. The bar graphs in inset compare the mean values of the reversal potential relative to RP in control females (black) and valproate females (red) at P20–23. A significant difference between control females and valproate females (∗∗p = 0.003) is reported.

in control mice. Therefore, the conductance shift is not impacted in this ASD model.

However, the developmental change in the GABAA receptor channel reversal potential was clearly more complicated in valproate-treated male mice (**Figure 4D**, black symbols) than in control male mice (red sigmoid, from **Figure 3C**). At P6–8 the reversal potential is significantly (p < 0.001, see graph in inset) more positive relative to RP (+20.6 ± 6.5 mV, n = 10) in valproate males compared to control males (+7.4 ± 4.6 mV, n = 9), decreases to a plateau value around +12 mV at P10, then reaches a negative value around P20 (sigmoid curve, X0 = P17.7 days). Thus, the switch from depolarizing to hyperpolarizing effect of GABA is delayed by 3 days in valproate male mice compared to control males. Therefore during a relatively long period, GABAergic currents are more depolarizing in valproate mice than in age matched control mice. Indeed, at P15–19 the mean reversal potential of the GABAA channel is still positive in valproate mice (+6.6 ± 4.4 mV, n = 17 compared to −15.9 ± 7.4 mV, n = 11 in control male mice, p < 0.001, see graph in inset). These differences are less evident in females (**Figure 4E**). Yet at P7 in valproate female mice the reversal potential of GABAA channels is around +20 mV relative to RP whereas it is about +5 mV relative to RP in control females. Furthermore, a comparison between the sigmoidal fit of the data obtained from valproate females (black curve, X0 22.7 days) and the sigmoidal fit of the data obtained in control females (red curve, from **Figure 3C**) reveals a delay of the GABA switch of almost 4 days. As a consequence, at P20-P23 the mean reversal potential of the GABAA channels is positive relative to RP (+5.8 ± 6.5 mV, n = 8) in valproate treated females but negative (−7.0 ± 3.0 mV, n = 8, p = 0.003) in control females (see inset **Figure 4D**).

FIGURE 5 | Distinct evolution of the GABA switch between conductance levels 1 and 2 in males, females and valproate-treated males. (A1) Dominant GABAA channel conductance from level 1 to level 2 as a function of the post–natal age in males (each point represents a single cell). Data are best fit with a sigmoidal hill function with four parameters (black line, a = 6.2, b = −1.33, X0 = 14.7 and Y0 = 11.1). (A2) Dominant GABAA channel conductance from level 1 to level 2 as a function of the post–natal age in females (each point represents a single cell). Data are best fit with a sigmoidal hill function with four parameters (red line, a = 7.8, b = −3.3, X0 = 17.7 and Y0 = 10.5). Note the presence of three channels with a conductance around 25 pS. (B1) Histogram distribution of the slope conductance of level 2 channels (pooled data from Figures 1B, 2B, 4A). The distribution is normal and fit by a Gaussian (blue line) with a peak value at 10.9 pA. (B2) Histogram distribution of the slope conductance of level 1 channels (pooled data from Figures 1B, 2B, 4A). The distribution is normal and fit by a Gaussian function (red line) with a peak value at 17.8 pA. (B3) Histogram distribution of the slope conductance of level 1 and level 2 channels (pooled data from B1 and B2). To determine the limit of conductance between level 2 (blue line) and level 1 (red line) the Gaussian fits obtained in (A1, A2) were superimposed. The limit between levels 1 and 2 is 14 pS. Channels with conductances ≥14 pS are classified as level 1 channels and channels with conductances <14 pS are classified as level 2 channels. (C1) Data from male mice. The bar plots show the % of level 1 channels (conductance ≥14 pS) at increasing post-natal age ranges, as indicated on top of each bar. (C2) Data from female mice. The bar plots show the % of level 1 channels (conductance ≥14 pS) at increasing post- natal age ranges (the same range as in C1). (C3) Percentage of conductance level 1 as the function of the age in male mice (black symbols) female mice (red symbols) and valproate treated male mice (blue symbols). Data for control males (black curve) and females (red curve) are best fit using a sigmoidal function with the following parameters: a = 99, b = −1.9, X0 = 15.3 and Y0 = 3.2 for males; and a = 100, b = −2.42, X0 = 18.5, Y0 = 0.43 for females.

# Sex Dependence of the GABA Conductance Shift: Parallel Development of GABAA Channel Reversal Potential and Conductance

To determine whether the GABAA channel conductance and the reversal potential switch occur in parallel, we compared the development of the GABAA channels conductance in males and females (**Figure 5**). As a first approach, the slope conductance of the dominant GABAA channel was plotted as a function of age in males (**Figure 5A1**) and females (**Figure 5A2**). The evolution of the conductance could be fitted in both sexes by a sigmoidal function (black and red lines) but with distinct parameters (slope = −1.33, X0 = 14.7 days for males and slope = −3.3, X0 = 17.7 days for females) indicating that the shift in conductance is slower and later in females compared

to males. In both sexes, a single channel with larger slope conductance (25–27 pS, not included in the fit) was recorded in three patches and may indicate the presence of a rarelyoccurring third conductance level of GABAA receptor channel, as previously described for cerebellar granule cells (Brickley et al., 1999).

trace, control conditions (in the presence of NBQX), lower trace in the presence of Isoguvacine. As in the male, isoguvacine increased firing

frequency and progressively reduced spike amplitude.

In order to quantitatively analyze the conductance switch, we determined the conductance ranges of the two levels. We first constructed separate histogram distributions of conductance levels 1 and 2 channels by pooling conductance values of GABAA channels recorded in both sexes between P5–12 for conductance level 1 (**Figure 5B1**) and conductance values of GABAA channels recorded at P28–45 in males and P26–36 in females for conductance level 2 (**Figure 5B2**). The distributions are fitted with distinct Gaussian functions, giving a peak of 17.8 pS for level 1 and 10.9 pS for level 2. Next, we combined the two distributions and associated fits (**Figure 5B3**), to show that the two fits intercept at a value of 13.9 pS. From this observation, we classified channels with a slope conductance of 14–22 pS as level 1, and channels with a slope conductance <14 pS as level 2.

Then, we determined the developmental alterations of level 1 channel percentage at different post-natal ages (**Figures 5C1–C3**). In both sexes the proportion of patches displaying dominant level 1 GABAA receptor channel decreases with age; but the histograms in **Figures 5C1,C2** show that this drop in level 1 conductance (defined as presence in less than 50% of patches) occurred earlier in males (P14; **Figure 5C1**) than in females (P17; **Figure 5C2**). Therefore, during a transitional period (around P15–20), PCs in females have a higher-conductance dominant channel than male PCs.

The evolution of the proportion of level 1 GABAA channel as a function of age (**Figure 5C3**) is similar in both males and females to the evolution of the chloride reversal potential. Both curves are sigmoidal (four parameters) with similar X0 s (14.4 days for the level 1 plot and 15.3 days for the chloride reversal potential plot in males; and 18.5 days for both the level 1 plot and the chloride reversal potential in females). Interestingly the switch of GABAA channel conductance is also altered in valproate male mice with 25% of level 1 channel dominant in young adult animals.

### Sexual Dimorphism of the GABA Excitatory/Inhibitory Shift

Finally, we determined possible links between the chloride driving force and the excitatory or inhibitory actions of the GABAA agonist isoguvacine on PC firing, with or without NBQX. Similar effects of isoguvacine were observed in the presence of NBQX suggesting that glutamatergic neurotransmission was not involved (not shown; n = 80). When the reversal potential of the GABAA channel was positive to RP (at P6), isoguvacine increased the firing rate in males (**Figure 6A**, n = 3 out of 5) and in females (**Figure 6B**, n = 3 out of 6). This effect was accompanied by a progressive decrease in the size of the spikes due to the sodium channel inactivation produced by a strong depolarization. In the remaining recordings, PCs were silent and isoguvacine applications led to a depolarization. In older PCs, when the reversal potential of the GABAA channels is negative to RP, isoguvacine decreased spontaneous activity of PCs, without altering the spike amplitude in both males (P27; **Figure 7A**) and females (P25; **Figure 7B**).

Additional experiments confirmed the sex and age difference in the actions of isoguvacine. In males at P15, after 2 min isoguvacine application, PC firing was blocked by 75 ± 35% in all cells tested (n = 4; example in **Figure 8A**), and completely blocked (100%) for PCs recorded from P16 to P29 (n = 14; example in **Figure 7A**). However in females at P17, 2 min of isoguvacine increased PC firing frequency by 220 ± 161%, accompanied by a reduction in spike amplitude in 5 out of 13 PCs (example in **Figure 8B**), sometimes preceded by a temporary reduction in the firing frequency (n = 2, including the example in **Figure 8B**). Isoguvacine had no effect on firing frequency in two PCs from female mice, and inhibited firing frequency in seven of these

about 7 Hz; isoguvacine abolishes the spiking activity of the cell; and this effect is reversible upon washout. Note that spike amplitude is not reduced before the silencing of the cell. (B) As in panel (A), but from a female P25 mouse. Upper panel: current traces recorded in control conditions, in the presence of Isoguvacine, and following washout. Isoguvacine silences the cell, and activity recovers after washout. Lower panel: data from the same cell, showing instantaneous action potential frequency (left) and action potential amplitude (right) over time. Isoguvacine is applied as indicated by the red line. The control spike frequency is around 30 Hz. In presence of isoguvacine the spiking activity of the cell is almost abolished, but the spike amplitude is not altered.

PCs. The excitatory effect of isoguvacine continued to be seen in female PCs (n = 2 out 5) until P19; then at P20–25 isoguvacine had only inhibitory effects (87% ± 8% inhibition, n = 6).

To summarize the effects of isoguvacine on spontaneous PC, firing we assigned a value of +1 when isoguvacine produced an increased firing frequency, 0.5 when the increase was preceded by a short inhibition of firing, 0 when isoguvacine had no effect and −1 when isoguvacine inhibited firing. **Figure 9** shows the mean valence as a function of age for PCs from male and female mice, and clearly illustrates the different time course of the shift to inhibitory isoguvacine effects between male and females.

# DISCUSSION

# Development of GABAA Receptor Channel Subunit Composition in Cerebellar PCs

Somatic GABAA receptor channels in adult vertebrate neurons have complex gating and multiple conductance states ranging from 7 to 36 pS, with a single dominant conductance state in each cell type (Newland et al., 1991; Robello et al., 1998; Brickley et al., 1999; Mortensen and Smart, 2006). We found three main conductance levels in PC, with conductances of 7–12 pS, 15–18 pS, and occasionally 25–28 pS. Similar conductance values are observed in adult cerebellar granule cells, but with a shift of the dominant conductance from medium (15–18 pS) to small (7–12 pS) during maturation. This shift likely reflects a reorganization of GABAA channel subunit composition (Moss et al., 1990; Verdoorn et al., 1990; Mortensen and Smart, 2006). Although it is difficult to equate GABA conductance levels with specific subunits combinations (Brickley et al., 1999), it is likely that in young-adult PCs, receptors lack the γ subunit (Verdoorn et al., 1990; Fisher and Macdonald, 1997; Amato et al., 1999; Hörtnagl et al., 2013). The dominant medium conductance level (18 pS) in newborn PCs suggests the presence of an α1β1 GABAA receptor at this developmental stage (Moss et al., 1990). The very rare occurrence of the 25–28 pS conducting level is may be related to a low level of expression of the γ2 subunit (Mortensen and Smart, 2006). Our observations suggest a subunit reorganization during development (but see

firing is reduced in the presence of isoguvacine and recovers after washout. Lower panel: data from the same cell showing instantaneous action potential frequency with time (left) and action potential amplitude with time (right). Isoguvacine application is indicated by the red line. Control firing frequency is about 30 Hz; isoguvacine abolishes firing, and this effect is reversible. Note that the abolition of the discharge is not linked to a decrease in spike amplitude. (B) As in panel (A), but from a P17 female mouse. Upper part: current traces recorded in control conditions, with isoguvacine, and after washout. Regular action potential firing is reduced in the presence of isoguvacine and recovers after washout. Lower graphs, from the same cell, show instantaneous action potential frequency with time (left) and action potential amplitude with time (right) Isoguvacine application is indicated by the red line. Control firing frequency is about 30 Hz; with isoguvacine application, a brief inhibition of the firing frequency is followed by an increase, up to 100 Hz; this effect is reversible. Isoguvacine also reversibly decreased spike amplitude.

Nadler et al., 1994; Haghir et al., 2013) but further experiments are required to determine if this is associated with changes in IPSCs kinetics.

#### Sexual Dimorphism of GABA Driving Force in Wild Type Mice and in the Valproate Model of Autism

In a wide range of brain structures, the polarity of GABAA receptor effects is regulated by the expression and activity of two chloride membrane transporters: NKCC1 and KCC2. NKCC1 is expressed in immature neurons and actively transports chloride ions into the cell; KCC2 is expressed in mature neurons and actively transports chloride ions out of the cell. This changes in transporter expression leads to depolarizing effects of GABA in immature neurons and hyperpolarizing effects in mature neurons (Ben-Ari, 2014; Watanabe and Fukuda, 2015; Raimondo et al., 2017) with a developmental switch around birth (Tyzio et al., 2014; Watanabe and Fukuda, 2015). However, PCs do not express NKCC1 (Mikawa et al., 2002), suggesting that alternative mechanisms are required to regulate chloride homeostasis, including chloride channels (Zhang et al., 2015; Rahmati et al., 2016) and KCC2. The developmental GABA switch occurs in PCs at the end of the second post-natal week (see also Eilers et al., 2001), when profound morphological and functional alterations occur in PCs (Dusart and Flamant, 2012). The timing of the switch is strongly dependent on the specific cerebral structure; it is similar in the antero-ventral cochlear nucleus (Song et al., 2012) and substantia nigra pars reticulata (SNR; Kyrozis et al., 2006), and earlier in other brain structures (Stein et al., 2004; Tyzio et al., 2007; Allain et al., 2011; Witte et al., 2014). Our results also suggest that the regulation of chloride homeostasis is sex dependent, since it is delayed in female mice. A similar sex dependence of chloride homeostasis has been observed in the substantia nigra (Galanopoulou, 2008).

The mechanisms underlying these sex differences are poorly understood. The testosterone surge in males that occurs during late gestation and shortly after birth (Dean and McCarthy, 2008) has been suggested to promote the expression and the activity of NKCC1 and to depress the synthesis and activity of KCC2 (Waddell and McCarthy, 2012), and thus make the GABA switch occurs later in males than in females. This effect has been shown in hippocampal cultures (Nuñez and McCarthy, 2007), SNR neurons in acute slices (Kyrozis et al., 2006), and embryonic hypothalamic neurons in culture (Mir et al., 2017). However, the GABA switch in PCs is delayed in females, rather than in males, despite this testosterone effect on KCC2, suggesting that alternative mechanisms must be present. PCs synthesize estradiol from cholesterol in the neonatal period and the enzymes involved in this synthesis are developmentally regulated with a different profile in males and in females (Dean and McCarthy, 2008). In PCs, estradiol promotes BDNF-mediated dendritic growth, spinogenesis and synaptogenesis during neonatal life (Tsutsui et al., 2011), and BDNF inhibits KCC2 activity (Huang et al., 2013). Furthermore, the development of the cerebellar expression of insulin receptors (IR) and insulin-like growth factor-1 receptors (IGF-1R) differs in males and females. IR and IGF-1R expression increase between P0–P7 and is down regulated at P14 in males, whereas in females IR and IGF-1R are stable between P0–P7 and up-regulated at P14 (Haghir et al., 2013). IGF-1 accelerates the developmental switch between NKCC1 and KCC2 chloride transporters in the visual cortex (Baroncelli et al., 2017). Collectively, these observations converge to suggest multiple mechanisms to explain the differences between chloride homeostasis regulation in males and females.

Using the valproate mouse model of autism, we observed that the chloride reversal potential in P5–6 valproate males is more positive than in controls (by about 10–20 mV), and that the GABA switch is delayed by 3–4 days in both males and females, extending this experimental model to the cerebellum. Increased excitatory-inhibitory ratio is suggested to underlie the pathogenesis of autism and successful clinical trials have shown that the NKCC1 antagonist bumetanide, which restores low (Cl−)I and GABAergic inhibition, also attenuates the severity of autism (Rubenstein and Merzenich, 2003; Lemonnier and Ben-Ari, 2010; Tyzio et al., 2014; Lee et al., 2016; Uzunova et al., 2016; Lemonnier et al., 2017). It is now accepted that the cerebellum is involved in higher order functions (Schmahmann and Sherman, 1998) including perceptual processes (Baumann et al., 2015), in addition to its roles in balance, posture, and motor control. The cerebellum has been implicated in many psychiatric disorders (Phillips et al., 2015) and cerebellar abnormality is associated with ASD (Tsai et al., 2012). Several reports in individuals with autism or animal models of autism describe alterations in PC density and properties (Skefos et al., 2014; Tsai, 2016). Thus, chloride gradient modifications and alterations in the GABA switch in cerebellar PCs that we show here may be particularly relevant for further studies of autism and the role of the cerebellum.

### Parallel Development of GABAA Receptor Channel Subunit Composition and Chloride Gradient

Our data show that the kinetics of the shift of the dominant conductance levels during development, from medium (15–18 pS) to small conductance (7–12 pS), are similar to the kinetics of the development of the chloride gradient. Both phenomena followed a similar sigmoidal process, are delayed in females, and are altered in valproate mice, suggesting interactions between the subunit composition of GABAA channels and intracellular chloride concentration (Succol et al., 2012).

Interestingly, in PCs, it has been suggested that an interaction between KCC2 and a specific subunit of the GABAA receptor represents a fundamental mechanism of regulation of GABAergic synapses (Huang et al., 2013). Furthermore, membrane expression of the β3 GABAA receptor subunit at different postnatal developmental stages is observed in rats exposed prenatally to valproic acid (Li et al., 2017). Finally, alterations in the efficacy of neuronal inhibition mediated by GABAA receptors containing β3 subunits have been implicated in autism (Vien et al., 2015). Collectively, these observations raise the possibility of parallel convergent developmental alterations of GABAergic signals, linking chloride gradient and subunit composition.

# CONCLUSIONS AND PERSPECTIVES

The conductance properties of GABAA receptors as well as the chloride equilibrium potential and the polarity of GABA effects in PCs are developmentally regulated in parallel, occurring later than in other brain structures, which is coherent with the delayed development of cerebellar PCs to the postnatal period. The switch is also delayed by 4 days in females, indicating sexual dimorphism in keeping with recently reported sex differences in synaptic excitation, inhibition and intrinsic properties (Mercer et al., 2016). Finally, we show that the valproate model of autism produced a shift in the chloride gradient in developing PCs from males and females, providing important potential implications for the cerebellum in the pathology of autism.

#### REFERENCES

#### AUTHOR CONTRIBUTIONS

SR, AL, YB-A, BP and J-LB wrote the manuscript. SR, AL and J-LB designed the research, performed the experiments and analyzed the data.

developing cerebellum. Neuron 38, 785–796. doi: 10.1016/s0896-6273(03) 00298-8

neurones of the rat. J. Physiol. 432, 203–233. doi: 10.1113/jphysiol.1991. sp018382

in the purkinje cell. Front. Endocrinol. 2:61. doi: 10.3389/fendo.2011. 00061

**Conflict of Interest Statement**: YB-A is CEO of Neurochlore, a biotech company devoted to treat autism.

The remaining 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 © 2018 Roux, Lohof, Ben-Ari, Poulain and Bossu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) 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.

# Intracellular Dynamics in Cuneate Nucleus Neurons Support Self-Stabilizing Learning of Generalizable Tactile Representations

Udaya B. Rongala1 , Anton Spanne2 , Alberto Mazzoni 1 , Fredrik Bengtsson2 , Calogero M. Oddo1 and Henrik Jörntell 2 \*

1 The BioRobotics Institute, Scuola Superiore Sant'Anna, Pisa, Italy, 2 Section for Neurobiology, Department of Experimental Medical Sciences, Biomedical Center, Lund University, Lund, Sweden

How the brain represents the external world is an unresolved issue for neuroscience, which could provide fundamental insights into brain circuitry operation and solutions for artificial intelligence and robotics. The neurons of the cuneate nucleus form the first interface for the sense of touch in the brain. They were previously shown to have a highly skewed synaptic weight distribution for tactile primary afferent inputs, suggesting that their connectivity is strongly shaped by learning. Here we first characterized the intracellular dynamics and inhibitory synaptic inputs of cuneate neurons in vivo and modeled their integration of tactile sensory inputs. We then replaced the tactile inputs with input from a sensorized bionic fingertip and modeled the learning-induced representations that emerged from varied sensory experiences. The model reproduced both the intrinsic membrane dynamics and the synaptic weight distributions observed in cuneate neurons in vivo. In terms of higher level model properties, individual cuneate neurons learnt to identify specific sets of correlated sensors, which at the population level resulted in a decomposition of the sensor space into its recurring high-dimensional components. Such vector components could be applied to identify both past and novel sensory experiences and likely correspond to the fundamental haptic input features these neurons encode in vivo. In addition, we show that the cuneate learning architecture is robust to a wide range of intrinsic parameter settings due to the neuronal intrinsic dynamics. Therefore, the architecture is a potentially generic solution for forming versatile representations of the external world in different sensor systems.

Keywords: cuneate nucleus, neurophysiology, neuronal plasticity, Intrinsic dynamics, tactile, touch, synaptic integration

# INTRODUCTION

The problem of how to represent a complex external world to support non-trivial versatility of action has a deadening presence both in neuroscience (Loeb and Fishel, 2014; Spanne and Jorntell, 2015), robotics and artificial intelligence (AI) systems (Service, 2014). For neuroscience, the issue is closely associated with the understanding of the brain—without

#### Edited by:

Philippe Isope, Centre National de la Recherche Scientifique (CNRS), France

#### Reviewed by:

Jesus Garrido, Universidad de Granada, Spain Natasha Alexandra Cayco Gajic, University College London, United Kingdom

#### \*Correspondence:

Henrik Jörntell<EMAIL_ADDRESS>
Received: 26 March 2018 Accepted: 26 June 2018 Published: 31 July 2018

#### Citation:

Rongala UB, Spanne A, Mazzoni A, Bengtsson F, Oddo CM and Jörntell H (2018) Intracellular Dynamics in Cuneate Nucleus Neurons Support Self-Stabilizing Learning of Generalizable Tactile Representations. Front. Cell. Neurosci. 12:210. doi: 10.3389/fncel.2018.00210

**19**

knowledge of how information of the world is represented in neuronal circuitry, it is difficult to decipher its functional mechanisms. An important related issue is how biological systems can generalize previously learnt sensory experiences to apply them to the interpretation of novel contexts—lack of generalization capability is a limitation in classical pattern recognition systems (Spanne and Jorntell, 2015), to which AI and DNN systems has an ancestry, and likely an important reason why such systems can be comparatively easily fooled (Nguyen et al., 2015). As versatility of interaction with the external world is a hallmark of brain function, an important question is how that versatility can be supported.

The skin is an interface that directly interacts with the physical world, using 10,000's of tactile sensors (Johansson and Flanagan, 2009). Current models of the organization of tactile neural systems to a large extent build on assumptions that the brain needs to identify edges, shapes or other physical parameters that human external observers may deem important to represent (Pei et al., 2011; Sathian, 2016). Using sparse-coding interpretations of neural coding combined with grandmother neuron-like theory (Barlow, 1972) such a system can be expected to work in the classification of a large range of tactile experiences. However, classifying systems based on pattern recognition can suffer from problems with generalization, i.e., where learning from one situation can be scaled or adapted to apply to new ones (Spanne and Jorntell, 2015). An alternative mode of representation of tactile information would be one that automatically arises out of experienced interactions through learning. Indeed, during early development, mammalians generate seemingly random movements and interactions with the environment, which would engage a wide set of sensors from large parts of the body (Shao et al., 2016) and play a crucial role for development (Forssberg et al., 1995; Petersson et al., 2003; Blumberg et al., 2013). Such interactions result in spatiotemporal patterns of skin sensor activation, that depend on, and therefore abstract, the properties of the objects we touch, the laws of physics describing the interactions that can take place (Hayward, 2011), the types of movement we make, the mechanical effects inducible in the skin and how the tactile sensors are tuned to them (Jörntell et al., 2014). The available set of sensors will respond across these conditions and their activations will have specific relationships depending on the condition. Hence, rather than viewing brain integration of tactile sensors as occurring in a pixel-by-pixel fashion, we here consider the often overlooked fact that individual neurons integrate information from several sensors. It follows that what is being learnt is likely to involve the relationships between the activations of these sensors. To learn such relationships is here hypothesized to be an important component of being able to form representations of the external world that is applicable or generalizable to novel experiences.

As tactile inputs are first processed in the cuneate nucleus, before they reach the cortex, it is likely that the basic constraints on the brain's representation of the tactile part of the physical world are formed here. In vivo whole cell recordings from these neurons indicate that their synaptic weights for tactile inputs are highly skewed, which indicates that they are highly learnt (Bengtsson et al., 2013). The important question why that learning occurs recently found a possible answer, when the cuneate neurons were found to code for unique combinations of the fundamental haptic input features (Jörntell et al., 2014), which tentatively correspond to the dimensions or the vector decompositions of the contact mechanics effects arising between two objects (Hayward, 2011). Here, we emulated the learning that would arise from of a set of varied sensory experiences given the biological constraints provided by the recording data we obtained on the intrinsic membrane dynamics and the synaptic inputs of the neurons of the cuneate nucleus. We find that the main effect of the cuneate learning is a utility-based decomposition of the tactile sensory space into vector components that made it possible to generalize the learning to novel tactile experiences. This is a different form of representation of sensory input data than a direct identification of sensory input patterns, which is a potential explanation for the large versatility in the identification of sensor input data in biological systems.

# METHODS

All biological data was obtained under acute conditions identical to those of a previous study on the cuneate nucleus in vivo (Bengtsson et al., 2013). Briefly, adult cats of both sexes were initially prepared under propofol anesthesia and subsequently decerebrated to allow cuneate neuronal recordings under nonanesthetized conditions. This study was carried out in accordance with the recommendations of Malmö-Lund Animal Research Ethics Committee. All experimental procedures were approved in advance by the Malmö/Lund Animal Research Ethics Committee (permit number and approval-ID: M32-09).

This section contains four main parts. First, we describe how recordings were made from projection neurons and inhibitory interneurons of the cuneate nucleus using the in vivo whole cell patch clamp technique. Secondly, the recorded characteristics of the cuneate projection neurons and the inhibition from the interneurons were approximated by constructing a model of individual cuneate projection neurons (CNs) and their afferent network. Thirdly, the responses of a population of skin sensors that provided synaptic input to the CN network across a range of different real world stimuli were generated by a bionic fingertip. Fourthly, the CN synaptic learning process for skin sensor input was inferred from generic neuronal learning mechanisms in vivo and our estimation of intracellular calcium responses. The theoretical basis for the construction of the CN model and its learning process is also provided for each step.

# Biological Data

Briefly, under initial deep propofol anesthesia, adult cats were prepared for acute in vivo recordings with artificial respiration, strict monitoring of blood pressure, end-expiratory carbon dioxide and temperature. Thereafter, the animals were decerebrated and the anesthesia discontinued. To monitor the level of anesthesia before decerebration, we continuously measured the blood pressure and verified the absence of withdrawal reflexes to noxious stimulation. To monitor the state of the preparation after the decerebration, we in addition made EEG recordings from intact parts of the neocortex. EEG recordings were characterized by a 1–4 Hz oscillatory activity that was periodically interrupted by large-amplitude 7–14 Hz spindle oscillations lasting for 0.5 s or more. Such EEG patterns are associated with deep stages of sleep (Niedermeyer and Da Silva, 2005). The EEG activity and the blood pressure remained stable, also on noxious stimulation, throughout the experiments. Mounting in a stereotaxic frame, drainage of cerebrospinal fluid, pneumothorax and clamping the spinal processes of a few cervical and lumbar vertebral bodies served to increase the mechanical stability of the preparation.

In vivo whole cell patch clamp recordings were made in the exposed cuneate nucleus. The recorded neurons were identified with respect to the location of their excitatory receptive field on the skin (Bengtsson et al., 2013). Stimulation in relation to the location of this receptive field was done using forcetime controlled touch as well as by local electrical activation of skin afferents. Current injections to characterize the intrinsic membrane responses were made through the recording pipette. All intracellular signals were analyzed off-line using custommade software. Identification of unitary IPSPs was made using a tailored template matching routine with manually constructed templates. The templates consisted of a time series of voltage ranges into which the signal had to fit to be counted. The template was adjusted to allow the detection of gradually smaller events until it started to include events that on visual inspection did not resemble the event time course of IPSPs evoked by skin stimulation (Bengtsson et al., 2013). Hence, the similarity of time course between evoked and unitary IPSPs was taken as an indication that they were derived from the same type of inhibitory synapses.

A difference with our previous analysis of EPSPs (Bengtsson et al., 2013), was that the voltage deflections of the unitary IPSPs were much smaller. Thereby, the noise of the recording system prevented the detection of IPSPs of amplitudes below a certain level. Therefore, our reported values of the median peak unitary IPSP amplitudes for each neuron are most likely overestimates. For population data, we report the mean and standard deviation of the median IPSP amplitudes recorded for at least 100 spontaneous IPSPs in each of the 15 recorded cuneate projection neurons.

#### Modeling Design

Based on our present biological observations, as well as previously published data on cuneate neurons (Bengtsson et al., 2013) and generic neuronal physiology and plasticity, we aimed to simulate how these parameters and processes could support learning when brought together in a functional system. Our simulation consisted of two main components,

(a) A bionic fingertip covered with silicon skin and equipped with tactile sensors that transduced local mechanical skin strain patterns into spatiotemporal patterns of sensor spike output data in response to physical interactions with external objects/surfaces. This was the counterpart of the tactile primary afferents (PAs) in biological systems. The bionic fingertip provided the important feature of a system of PA sensors where there is a consistent relationship between the activation of the different sensors across conditions or tactile experiences. As discussed in the introduction, rather than viewing brain integration of tactile sensors as occurring in a pixel-by-pixel fashion, we here consider the often overlooked fact that each cuneate nucleus neurons integrate information from several sensors. It follows that what is being learnt is likely to involve the relationships between the activations of these sensors. To learn such relationships is here hypothesized to be an important component of being able to form representations of the external world that is applicable to novel experiences.

(b) A simulated cuneate nucleus neuronal network containing a set of cuneate neurons receiving synaptic input from those PAs. Each model cuneate projection neuron (CN) was simulated individually, capturing the intrinsic responsiveness observed in vivo and the generic subcellular processes involved in synaptic plasticity. The accumulative effects of the synaptic plasticity during the learning process were evaluated by the changes in synaptic weights of the PA inputs in each CN.

In the account that follows, we present the design of the simulation starting with the low level neuronal properties, the dynamic model describing the intrinsic responsiveness of the cuneate neurons, synaptic weight initialization, the sensory input data and the bionic fingertip used to generate it, the learning process and the evaluation of the end effect of the learning process. The design and performance of the bionic fingertip has already been reported in previous papers, therefore we merely present the specific adaptations made for the present simulation. The modeling software is available on: https://figshare.com/ projects/Artificial\_Cuneate\_Neuron\_Model/34886.

#### Network Connectivity

Based on biological data of cuneate neurons (Bengtsson et al., 2013; **Figures 1**, **2**), our network model (**Figure 3A**) was comprised of sensory afferents projecting as excitatory synapses on each individual model CN. In addition to the excitatory synapses, and based on our biological observations (**Figure 1**), the sensory afferents also drove local inhibitory interneurons that provided inhibitory synapses to the CNs (**Figure 3A**). The model CN had equal numbers of excitatory and inhibitory synapses. The synaptic weights of these synapses were given initial starting weights (or seed weights, see below) that were gradually learnt on a trial-by-trial basis.

The simulated network model had 80 sensory input channels (PAs) that innervated each individual CN. This number of afferents was less than biological estimates, which suggest in the order of 1,000's of PAs per CN (Bengtsson et al., 2013). However, the lower number of simulated afferents is still realistic because out of these 1,000 synapses, many are likely to represent anatomically present synaptic inputs from PA axons mediating input from several different fingers, which in the adult animal mostly provide "silent" synaptic inputs (Weinberg et al., 1990; Bengtsson et al., 2013). In contrast, our inputs were generated from the tip of one finger alone. Synaptic inhibition was simulated as being provided by 80 independent interneurons that were each directly activated by one out of the 80 PA afferents available in our simulated system. Because all unitary inhibitory synaptic inputs were found to be of comparable, low weights in the biological cuneate neuron (**Figure 1**), inhibitory synapses were here simplified to one lump inhibitory synapse per CN (total inhibitory synaptic weight, **Figure 3**). Hence, interneurons were not explicitly simulated, but instead the spiking activity of each PA afferent was fed directly to an inhibitory synapse on the CN.

#### Neuron Model

The model cuneate projection neuron (CN) was implemented as a conductance based Exponential Integrate and Fire (EIF) model (Fourcaud-Trocmé et al., 2003), which can recreate the fast dynamics (∼1 ms timescale) of neuronal spike generation. In addition to the basic EIF model, voltage sensitive calcium channels and calcium dependent potassium channels were also modeled in order to recreate intermediate cuneate neuron dynamics (∼10 ms timescale) observed in vivo (**Figures 1E,F**). The complete dynamics of the CN membrane potential are given by:

$$C\_m \frac{dV\_m}{dt} = \left[I\_L + I\_{spike} + I\_{ion} + I\_{ext} + I\_{syn}\right] \tag{1}$$

where Cm isthe membrane capacitance, IL = −gL (Vm−EL)is the leak current, Ispike is the spike currents (fast dynamics), Iion is the ion channel currents (intermediate dynamics), Isyn is the synaptic input currents and Iext denotes external injected currents that were used to evaluate the intrinsic CN neuron responsiveness to current step commands (**Figure 1F**). The model values for these and other parameters are given in **Table 2**. The membrane resistance and time constant/capacitance were within the range of values recorded in vivo (Bengtsson et al., 2013), whereas other intrinsic parameter values were chosen through an optimization process (see below).

#### Neuron Model-Fast Dynamics

The spike current (Ispike) is generated using a basic EIF model (Equation 2) to achieve the fast dynamics and recreate the initiation of the action potential (Ostojic and Brunel, 2011).

$$I\_{spike} = \overline{\mathcal{g}}\_L \Delta\_T \exp\left(\frac{V\_m - V\_t}{\Delta\_T}\right) \tag{2}$$

As the depolarization reaches a threshold speed, an action potential is said to be generated, and the membrane potential is immediately reset by a fast afterhyperpolarization (AHP) to the membrane potential at which the spike was initiated, i.e., mimicking the fast spike AHP of the recorded cuneate projection neurons (see **Figure 2A**).

#### Neuron Model-Intermediate Dynamics

The intermediate dynamics include currents from additional ion channels (Iion, Equation 1) that are not directly involved in forming action potentials, but have more of a modulating role in episodes leading up to the generation of action potential and the episodes between action potentials when the synaptic input activity is high. The intermediate dynamics of the model were optimized to mimic the reactive conductances that could generate the types of responses to current injections we recorded in the cuneate neurons in vivo (**Figures 1E,F**). Such responses, i.e., post-inhibitory rebound responses and a tendency to generate bursts of action potential firing, have been observed in other neuron types (Llinás and Jahnsen, 1982; Huguenard, 1996; Molineux et al., 2008) and has at least partly been attributed to by low-threshold voltage gated calcium channels (LVA) and calcium-activated potassium channels (CAP). Hence, Iion can be divided according to:

$$I\_{\rm ion} = \; I\_{\rm Ca} + \; I\_{\rm K} \tag{3}$$

where ICa is the current through LVA channels or equivalent channels and IK the current through the CAP channels. These are modeled as two separate pools of ion channels (Saarinen et al., 2008) according to:

$$I\_{Ca} = -\overline{\mathfrak{g}}\_{Ca} \mathfrak{x}\_{Ca,a}^{3} \mathfrak{x}\_{Ca,i} \begin{pmatrix} V\_m - E\_{Ca} \end{pmatrix}$$

$$I\_K = -\overline{\mathfrak{g}}\_{K} \mathfrak{x}\_{K\_{Ca}}^{4} \mathfrak{x}\_{K\_{Vm}}^{4} (V\_m - E\_K) \tag{4}$$

where gCa and gK are the maximum conductances of the respective channels, ECa and EK are the reversal potentials of the respective ions, and xKCa , xKVm , xCa,a, xCa,i are the activity states of the channels (Saarinen et al., 2008).

The activation states of the LVA channels were modeled using the following differential equations:

$$\frac{d(\mathbf{x}\_{\text{C}\_{a,i}})}{dt} = \left(\overline{\mathbf{x}}\_{\text{C}\_{a,a}} \left(V\_m\right) - \mathbf{x}\_{\text{C}\_{a,i}}\right) / \tau\_{\text{C}\_{a,i}}$$

$$\frac{d(\mathbf{x}\_{\text{C}\_{a,i}})}{dt} = \left(\overline{\mathbf{x}}\_{\text{C}\_{a,i}} \left(V\_m\right) - \mathbf{x}\_{\text{C}\_{a,i}}\right) / \tau\_{\text{C}\_{a,i}}\tag{5}$$

where tau is the time constant at which the states move toward the voltage dependent equilibrium described by xCa (Vm). These equilibrium functions are two parameters sigmoid of the form:

$$\overline{\chi}\_{\text{Ca},a} \left( V\_m \right) = \left( 1 + \exp \frac{p\_{\text{Ca},a,1} - V\_m}{p\_{\text{Ca},a,2}} \right)^{-1}$$

$$\overline{\chi}\_{\text{Ca},i} \left( V\_m \right) = 1 - \left( 1 + \exp \frac{p\_{\text{Ca},i,1} - V\_m}{p\_{\text{Ca},i,2}} \right)^{-1} \tag{6}$$

Since the calcium sensitive subunits of the CAP channels are located on the inner surface of the cell membrane, the intracellular concentration of calcium (- Ca2+ ) is modeled for small volumes rather than as an overall concentration within the cell. The calcium concentration within the cell will change both due to ion channels through which calcium ions enter to the inside of the membrane, and diffusion of ions into the remaining intracellular volume. From this line of reasoning, and in accordance with Saarinen et al. (2008), Equation 7 is constructed

$$\begin{split} \frac{d(\left[Ca^{2+}\right])}{dt} &= B\_{Ca}\overline{\mathfrak{z}}C\_{a}\mathfrak{x}\_{\text{Ca}\_{a}}^{3}\mathfrak{x}\_{\text{Ca}\_{a}}\left(V\_{m} - E\_{\text{Ca}}\right) \\ &+ \left(\left[Ca^{2+}\right]\_{\text{rest}} - \left[Ca^{2+}\right]\right)/\mathfrak{r}\_{\left[Ca^{2+}\right]} \end{split} \tag{7}$$

The activation states of the CAP channels were modeled using the following differential equations:

$$\frac{d(\mathbf{x}\_{K\_{\rm Ca}})}{dt} = \left(\overline{\mathbf{x}}\_{K\_{\rm Ca}} \left(\left[Ca^{2+}\right]\right) - \left[\mathbf{x}\_{K\_{\rm Ca}}\right]/\tau\_{K\_{\rm Ca}}\right)$$

$$\frac{d(\mathbf{x}\_{K\_{\rm Mn}})}{dt} = \left(\overline{\mathbf{x}}\_{K\_{\rm Mn}} \left(V\_m\right) - \left[\mathbf{x}\_{K\_{\rm Mn}}\right]\right)/\tau\_{K\_{\rm Mn}}\tag{8}$$

where the two time constants (τKCa and τKVm ) indicate the times at which the states move toward the voltage dependent equilibrium described by xKVm and the calcium dependent equilibrium described by xKCa . These equilibrium functions are two parameters sigmoid of the form:

$$\overline{\chi}\_{K\_{Vm}}(V\_m) = \left(1 + \exp\frac{p\_{K\_{Vm,1}} - V\_m}{p\_{K\_{Vm,2}}}\right)^{-1}$$

$$\overline{\chi}\_{K\_{Ca}}\left([Ca^{2+}]\right) = 1 - \left(1 + \exp\frac{p\_{K\_{Ca,1}} - [Ca^{2+}]}{p\_{K\_{Ca,2}}}\right)^{-1} \tag{9}$$

The values for these parameters were chosen through an optimization procedure (see below) and are indicated in **Table 2**. Note that we observed a range of variance between cuneate neurons in vivo in terms of their intermediate dynamics. Our aim was to provide a simple model that at the same time could capture the main principles of the rebound and burst responses that we could demonstrate in our recorded cuneate neurons. Therefore, we do not expect that our parameter values have a direct correspondence with biophysical measures, and we do not expect to precisely capture the properties of the intermediate dynamic response of any single neuron (which likely would have required a larger set of parameters).

#### Neuron Model-Synaptic Inputs

The synaptic current (Isyn) through the cell membrane is the summated synaptic currents of the activated synapses. Each individual synapse (i) is activated by a primary afferent spike generated by a single sensor of the bionic fingertip. Once activated, this spike gives rise to a stereotyped time course of conductance injection at the synapse (Isyn) which is described by

$$I\_{\rm syn} = \text{g\_{\rm max}} \sum\_{i} w\_{\rm exc,i} \exp\left(-\tau \left(t - t^{\
u}\right)\right) \left(E\_{\rm rev,exc} - V\_m\right)$$

$$+ \text{g\_{\rm max}} \,\omega\_{\rm inh} \sum\_{i} \exp\left(-\tau \left(t - t^{\
u}\right)\right) \left(E\_{\rm rev,inh} - V\_m\right) \tag{10}$$

where Erev is the reversal potential of the type of synapse (Erev,exc or Erev,inh depending on whether the synapse is excitatory or inhibitory, see **Table 2**), Vm is the membrane potential and t ∗ is the time of activation of the synapse. Each spike in each sensory afferent was converted into a synaptic conductance in the simulated neuron. For each synapse, the peak amplitude of the synaptic response was determined by the product of their individual weight (wexc or winh) and the overall maximum synaptic conductance constant (gmax, see **Table 1**). Through gmax the relative leak conductance (i.e., the ratio of the synaptic and the leak conductances) could be adjusted to simulate cuneate TABLE 1 | Specification of the intrinsic CN parameter values used.

The values of the basic CN configuration, or default values, are indicated in bold.

neurons with different sizes. The time constant of the decay of the synaptic membrane potential responses, τ , was 6.4 ms for both excitatory and inhibitory synapses, in accordance with the time courses recorded in the cuneate nucleus neurons in vivo (Bengtsson et al., 2013; **Figure 1**). Note that as all the PA synapses of our system stayed well below 200 Hz of firing activity, we did not simulate any rate adaptation of the PA synapses as such adaptation in vivo primarily occurs at intervals shorter than 5 ms (Jörntell and Ekerot, 2006).

#### Optimization of the Neuronal Calcium Dynamic Model Against Measured Data

The complete model was optimized during three steps. During each simulation, the model was fed with six 100 ms current steps with amplitudes of 100, 200, 300, 400, 600, and 800 pA. The results from the six trials were then optimized against intracellular recordings where the cuneate neurons were fed with the same currents in vivo (**Figure 1F**).

The first step is to manually choose suitable initial parameters, using both previously known values to some of the parameters, and estimating others using trial and error simulations. The second step is to use the Nelder-Mead algorithm, but with an objective function where the simulated traces is compared to the measured traces. The third step also uses the Nelder-Mead algorithm, but with an objective function that measure the discrepancy between the action-potential timing in the simulated trace and the measured trace. As there is no guarantee that the simulated trace contains the same number of action potentials as the measured trace, discontinuities appear when the number of action potentials in the simulated trace change. The use of the Nelder–Mead method is motivated by that this is a commonly applied numerical method used to find the minimum or maximum of an objective function in a multidimensional space, in particular for nonlinear optimization problems. Equation 11 contains the complete objective function:

$$e = \left\{ \begin{array}{c} e\_t \rho\_{\mathcal{S}} \\ \end{array} \sum\_{T, \hat{T} \in \mathbb{T}} \left( \left| \hat{T} \right| - \left| T \right| \right)^2 & \exists T, \hat{T} \in \mathbb{T} \quad \left| \hat{T} \right| \neq \left| T \right| \\\\ e\_t & \forall T, \hat{T} \in \mathbb{T} \quad \left| \hat{T} \right| = \left| T \right| \end{array} \right. \tag{11}$$

#### TABLE 2 | Neuron model parameters.

where e is the objective function value, ωs a weight used to punish any discrepancy in the number of action potentials, and et the total time error between simulated and measured action potentials:

$$e\_t = \sum\_{T, \hat{T} \in \mathbb{T}} \left\{ \sum\_{\substack{t \in T \\ t \in \hat{T}}} \min\_{\hat{t} \in \hat{T}} \left| t - \hat{t} \right|^2 \quad \left| T \right| \ge \left| \hat{T} \right| \right. \tag{12}$$

where T is the set of all pairs of T, Tˆ, where T and Tˆ are the sets of all measured and simulated action potentials, respectively, during a single current step. The list of all optimized parameters used in the model is shown in **Table 2**.

#### Subsynaptic Local Calcium Activity

In the learning process, excitatory synaptic weight learning was driven by the calcium activity in the main compartment of the cuneate neuron (i.e., as calculated by the calcium dynamic model above) in combination with the calcium activity in the individual synapses. An essential component of this combination is the intensity of activation of the individual synapses. According to the learning rule that we used (see below), a synapse that fires at high frequency with a high degree of correlation with the main compartment total calcium activity (A Ca2+ tot = k[Ca2+], where k is an arbitrary constant that is here assumed to be 1), will be "rewarded" as due to the strong correlation with the learning signal A Ca2+ tot . Conversely, strong firing in a synapse in relation to low or zero A Ca2+ tot will be "punished" (i.e., similar to the classical BCM rule for Hebbian plasticity) (Bienenstock et al., 1982). Therefore, the local calcium time constants (τCa2+ loc , **Table 1**), defining the temporal properties of the calcium signal in the local space underneath each individual synapse, play a major role in the learning process (the local postsynaptic calcium activity can be considered an analogy with the calcium activity in a local dendritic spine; Koester and Sakmann, 1998; Tigaret et al., 2016). The learning rule critically depends on this time constant. For instance, if τCa2+ loc time constants are too high, the rewarding effects on synapses that have a high degree of correlation with the A Ca2+ loc will be lost. However, as there is no data on the relevant time constants in the cuneate neurons in vivo, we had to make assumptions of the values of this time constant. In order to avoid pitfalls in relation to this assumption, we studied a range of time constants for A Ca2+ loc (**Table 1**) during the CN learning process.

For each synapse, each input spike at time t ∗ contributes to the subsynaptic spine calcium concentration, the time course of which is given by the kernel (Mazzoni et al., 2008) of Equation 13,

$$\begin{aligned} \alpha &= \frac{\mathfrak{r}\_1}{\mathfrak{r}\_d - \mathfrak{r}\_r} \\ \Delta A\_{loc}^{\triangle \hat{\mathfrak{a}}^{2+}}(t) &= \alpha \ast \left[ \exp \left( -\frac{t - \mathfrak{r}\_l - t^\*}{\mathfrak{r}\_d} \right) - \exp \left( -\frac{t - \mathfrak{r}\_l - t^\*}{\mathfrak{r}\_r} \right) \right] \end{aligned} \tag{13}$$

In its basic configuration, the parameters describing the relative local calcium concentration (or activity), are the decay time τd = 12.5 ms and the rise time τr = 4 ms multiplied with a constant τ1 = 21 ms (which is a constant to calculate the ratio) for τCa2+ loc = 100% (**Table 1**). τl is the latency time which is zero in our case. The initial values chosen were derived from our assumption that the time course of the slow afterhyperpolarization of the cuneate neuron spike (**Figure 1B**), which is known to reflect the activation of calcium-dependent potassium channels, matched the time course of the calcium concentration induced in the synapse. This resulted in a somewhat faster but still comparable time course of spine calcium than reported for single spines in vitro (Tigaret et al., 2016), but the rise time of our subsynaptic calcium signal was slower than in previous simulations of the

properties of calcium-dependent learning (Graupner and Brunel, 2012). As the temporal properties of this calcium signal clearly were assumptions with large uncertainty, we tested a wide range of different values for these parameters (**Table 1**). In order to achieve supralinearity in the local calcium activity (Wang et al., 2000), we used an approximative approach of subtracting an offset from the local calcium signal (an offset corresponding to 75% of the peak activity of the single pulse activation was subtracted) and this was the value of the local calcium activity (A Ca2+ loc ). With this approach, repetitive activation of the same synapse resulted in a supralinear increase in the intensity of the local calcium activity depending on the frequency of the PA afferent activity for that synapse.

#### Synaptic Weights

The synaptic weight in our model was a gain factor ranging between 0.001 and 1, where 1 signified the maximum weight of a synapse. As there is no information on perinatal synaptic weights for the CNs, we needed to make assumptions. The first assumption was that synapses have different initial synaptic weights, or seed weights. The second assumption was that all synapses had initial non-zero up to medium strong synaptic seed weights. In our model, the distribution of synaptic weights across the primary afferent synapses of the CN model, or the initial excitatory synaptic weights (winit,exc), were normal distributions ranging between 0.001 and 0.5 across the 80 primary afferent inputs. We used 5 different pseudo-randomized initial weight distributions, referred to as "Seed weights 1– 5." Pseudo-randomization here implies that the distribution is randomized but that the same distribution is kept constant across different learning process runs. This had the advantage that the effects of specific intrinsic CN configurations (see "Variations of initial intrinsic CN parameters" below) could be tested for the same initial weight distributions (**Figure 12**).

Synaptic inhibition was simulated as being provided by 80 independent interneurons that were each directly activated by one out of the 80 PA afferents available in our simulated system. Each PA synapse had the same weight on its targeted inhibitory interneuron. The collective, or total, inhibitory synaptic weight (winit,inh) was initially set to 0.125 evenly distributed across all the inhibitory synapses between the interneurons and the CN neuron they contacted, meaning that each PA synapse provided equally weighted inhibition to the CN neuron (as suggested by our biological data accounted for in the Results).

#### Sensory Inputs

Our aim was to simulate the learning process in the cuneate neurons driven by simulated sensory activation of the PAs. The main idea was to simulate a varied set of interactions between the skin and the physical world, and let the resulting spatiotemporal patterns of skin sensor activation determine the outcome of the learning process in the cuneate neurons. Rather than designing an arbitrary set of spatiotemporal patterns of skin sensor activations, we wanted to use a physical model of a fingertip to generate the spatiotemporal spike patterns of PA input that the cuneate neurons of our model (CNs) learned from. This is because across the different kinds of interactions that the skin may experience, there may be relationships between the skin sensors that are not easily calculated across all conceivable conditions/interactions.

### Bionic Fingertip

To generate the sensor patterns, we used a set of touch protocols based on different stimulus conditions (Table S1). The sensory fingertip comprises a 2 × 2 array of Micro Electro Mechanical System (MEMS) based bio-inspired tactile sensors (Oddo et al., 2011) to generate spatio-temporal sensory input patterns (Figure S2). Each individual sensor comprises a four transducing piezoresistors mesh (totaling 16 sensory channels for each fingertip), arranged in cross-shaped fashion able to generate precise response magnitudes for both normal and tangential forces, applied across the surface of the sensor (Beccai et al., 2005; Oddo et al., 2007). However, only four of the 16 sensor channels were sufficiently dynamically sensitive to the range of forces arising for the stimuli used to be used in the present experiment. We created 80 PA sensory input channels from these four sensors' analog data (**Figure 4**) by multiplexing them with multiple neuron models and signal processing as explained below. Each of these 80 PA sensor spike output patterns were considered a unique PA signal provided as synaptic input to the CNs, where sensor#1 was used to create PA input #1-20, sensor #2 to create PA input #21-40 and so on.

The normalized analog output of the tactile sensor is fed as virtual current input to a fine-tuned Izhikevich neuron models (Izhikevich, 2003) in-order to achieve the spatiotemporal spike output (Figure S2) as described previously (Oddo et al., 2016; Rongala et al., 2017). The Izhikevich model was chosen in order to reproduce the adaptation dynamics, that is a characteristic of mechanoreceptors (Johansson and Flanagan, 2009). Per the Izhikevich model, the membrane potential v and the adaptation variable u were updated via the following nonlinear differential equations discretized using Euler's method:

$$\dot{\nu} = A\nu^2 + B\nu + C - \mu + \frac{I\_{input}}{C\_m} \tag{14}$$

$$\dot{u} = a(b\nu - u)$$

When the membrane potential reached the spike depolarization threshold of 30 mV, one spike was produced followed by a reset:

$$\text{if } \nu \ge 30 \, mV, \text{ then } \begin{cases} \quad \nu \gets \mathcal{c} \\ u \gets u + d \end{cases} \tag{15}$$

The A, B, C and the spiking threshold are the standard ones of the Izhikevich artificial neuron model, whereas the parameters a, b, c and d were selected (Table S2) to allow a regular spiking behavior.

From the existing 16 sensory channels, we consider 4 active channels (2, 5, 12, and 15) (Figure S2). We further derivated the analog sensor data to mimic both Merkel Slowly Adapting (SA) type I mechanoreceptors, sensitive to the sustained level of mechanical interactions, and Meissner Fast Adapting (FA) type I mechanoreceptors, mostly sensitive to dynamic changes (Vallbo and Johansson, 1984; Abraira and Ginty, 2013). For each physical stimulus, the stimulation was repeated five times and the corresponding consecutive five analog signals were considered as the response of a unique sensor. As we had four physical sensors, this procedure thus generated 20 analog unique sensory signals per physical stimulus. In addition, each analog sensor signal was derivated to obtain a total of 40 different analog sensory signals (**Figure 4**). Moreover, by implementing the Izhikevich neuron model (Izhikevich, 2003) with two different settings, we obtained a total of 80 different PA spike output patterns for each physical stimulus.

For the learning process, we used the sensor outputs obtained using five different physical stimuli (**Figure 4**), and in addition the PA spike output responses to five other physical stimuli were provided as "non-training" stimuli (see Results text).

#### Learning Protocol

In order to induce learning in the CNs, we used a series of 1,500 stimulus presentations. Each of the 1,500 stimuli hence corresponded to a spatiotemporal pattern of PA spike inputs (**Figure 4**), where all of the 80 PAs were activated during each stimulus presentation. The 1500 presentations were generated from the 5 physical stimulus presentations. Each stimulus presentation was hence repeated 300 times. Rather than feeding the system with identical spatiotemporal spike patterns for each of these 300 repetitions, we added a spike time jitter to all tactile sensor spike output responses. Gaussian noise with zero mean and standard deviation σ (σ = 5 ms) was added to the individual spike times of each PA. The motivation for this injection of noise rather than repeating identical patterns of PA spike responses was because we observed it to make the learning more robust. It also allowed an identification of the spread of learning outcome (end weight "landscapes") for each model setting (**Figures 5B,C**).

#### Excitatory Synaptic Weight Learning

During the learning process, the individual excitatory synaptic weights were gradually learnt, i.e., there was an alteration in the weight of each synapse. All the synaptic weights were updated after each stimulus presentation. The weight change for the individual excitatory synapse (i) is given by Equation 16.

$$\text{cov}\_{\text{exc},i} = \int\_{t\_0}^{t\_{\text{max}}} \left\{ \left( A\_{\text{tot}}^{\text{Ca}^{2+}} \left( t \right) - \left( A \text{vg}\_{A\_{\text{tot}}^{\text{Ca}^{2+}}} . Sym\_{EQ} \right) \right) \cdot A\_{\text{loc}}^{\text{Ca}^{2+}} \left( t \right) \right\} \, K \, . dt \tag{16}$$

The driving force for the net synaptic weight change (wcexc,i) is given by the integral of the correlation between main compartment total calcium (A Ca2+ tot ) and local calcium activity for each synapse (A Ca2+ loc ) from t0 to tmax (corresponding to the start and end of each stimulus presentation). The A Ca2+ tot is offset to zero by the Learning threshold (**Figure 3D**), where A Ca2+ tot above the Learning threshold is in the positive zone (potentiation) and below is the negative zone (depression) (similar to Graupner and Brunel, 2012). In other words, the main compartment calcium activity is a gate for the plastic reaction in the local synaptic compartment, deciding whether potentiation or depression should occur. Therefore, the product of the counterbalanced A Ca2+ tot and individual synapse A Ca2+ loc (**Figures 3E–G**) defines the net learning drive for each respective synapse as a function over time (**Figures 3F–H**). The value of the time integral of the net learning drive attained at tmax decides to what extent each specific synapse should be potentiated (**Figure 3F**) or depressed (**Figure 3H**). The strength of the potentiation/depression for that individual synapse is further multiplied with a constant given by a sigmoid function (**Figure 3I**) of the current synaptic weight (wexc,i), a constant which we called the synaptic weight compensation. This is motivated by that the insertion of a synaptic receptor ion channel in the synaptic membrane can be regarded as a chemical reaction, where the addition of a new channel will become less likely the higher the number of ion channels already inserted and vice versa. The sigmoid function is defined by S (t) = 1 1+e−t , with an arbitrarily chosen steepness gain of 0.005.

The Learning threshold is given by the product of the average total calcium (AvgA Ca2+ tot ) and synaptic equilibrium (SynEQ) where AvgA Ca2+ tot is the mean of A Ca2+ tot across the last three consecutive stimulus presentations. The averaging was required to avoid instability in the learning and could in principle correspond to the dampening effect of the time required in biological systems for protein synthesis and/or ion channel insertion. The SynEQ is a gain factor, used to attain a homeostatic synaptic plasticity (Turrigiano and Nelson, 2004; Turrigiano, 2011) to keep the sum of the excitatory synaptic weights in control. It is defined as a linear function of the total excitatory synaptic weight (P wexc), with a dual slope having zero point preset to SetPointP wexc (in our model, that point is set to 10). The slopes intercepts for this linear function are defined such that, when total excitatory synaptic weight is 10 the SynEQ is 1. When P wexc is >10, the slope is 0.12, whereas for P wexc < 10, the slope is 0.04. The differences in slopes were necessary to prevent neurons from being stuck in a depressed state, and to prevent unstable, rapid potentiation once above the set point. In principle, the slope differences can be regarded as corresponding to that the insertion and the removal of synaptic conductances in biology is handled by partly separate chemical reactions (Jörntell and Hansel, 2006).

#### Inhibitory Synaptic Weight Learning

In model CNs without adaptation of synaptic inhibition, we observed that synaptic weight changes in the PA synapses tended to make the CNs unstable, which degraded learning stability (**Figure 10**). Instability was attained gradually during the learning process and typically resulted in that either the CN had excessive calcium activity or no calcium activity at all. In such cases, the CN would often swing between these two extremities several times during the learning process. We found that it was important that the model system prevented the CNs from entering such unstable states and a main effective regulator was found to be gradual adaptation of the total inhibitory synaptic weight around a set point of the total calcium activity of the neuron.

Therefore, during the learning process, after each stimulus presentation, also the total inhibitory weight was adapted toward a set point, which was defined by the activity of the calcium spike rate, closely related to A Ca2+ tot . A high calcium activity during the stimulus presentations was countered by an increase in total inhibitory synaptic weight, whereas a low calcium activity was countered by a decrease in this weight. A plastic regulation of the weight of inhibitory synapses against postsynaptic calcium activity or postsynaptic calcium-dependent bursting activity has previously been found in various neurons across the brain (Kawaguchi et al., 2011; Hulme and Connelly, 2014; Lourenço et al., 2014). In our model, we found this feature to prevent overexcitation in the simulated CNs as the learning of the synaptic weights progressed. A function with dual linear slope was used to define the rate of inhibitory weight change with respect to the calcium spike rate. Inhibitory weight change is summed with the existing inhibitory synaptic weights (winh). The calcium activity set point (SetPointCa2+ ) (**Table 1**) is the critical factor that defines the calcium spike rate, around which the inhibitory weight adapts to keep the system in balance. These slope functions are built zeroing from the preset point (SetPointCa2+ ), where the range of weight change is set between −0.01 and 0.01 for calcium spike rates between 0 and 200 Hz, respectively. In order to dampen the rate of adaptation, we used a moving average of the calcium spike rate across the last three stimulus presentations as the input to the adaptation.

Adaptation of the total inhibitory synaptic weights was implemented because the change in the weights of the excitatory PA synapses on the CN neuron could lead to large changes in the levels of the LVA (voltage-activated calcium channel) activity, which in turn led to unstable learning (as in **Figure 10**). Adaptation of the inhibitory synaptic weights balanced the LVA activity and resulted in a maintained stable learning process. All adaptation of the weight of the synaptic inhibition was evenly distributed across all the inhibitory synapses. Calcium-dependent potentiation of inhibitory synapses has previously been found to be confined to the inhibitory synapses that were active during the calcium activity (Kawaguchi et al., 2011; Hulme and Connelly, 2014; Lourenço et al., 2014)–indeed, in our simulated set of inhibitory synapses, all of them were active under all stimulus presentations as described above.

#### Variations of Initial Intrinsic CN Parameters

In addition to testing the effects of excitatory synaptic seed weight distributions on the outcome of the learning process (see below) we also wanted to test the effects of intrinsic CN parameters. Moreover, the values of these parameters to some extent relied on assumptions, and we wanted to know the sensitivity of the model to these assumptions. Therefore, in separate simulation runs, some of the intrinsic CN parameters were varied as described by **Table 1**. All of the initial parameters were kept constant during the learning process.

# MODEL DATA ANALYSIS

#### Correlation Index

The correlation index measure (**Figure 5D**) was used to quantify the degree of correlation between the spiking activities of two or more PAs while being provided with sensory inputs. We compared PAs whose synapses after learning ended up as High End Weight (HEW) synapses with high seed weight synapses that ended up as Low End Weight (LEW) synapses.

To compute this measure, we considered 10 HEW (green bars in **Figure 5C**) and 10 LEW synapses (red bars in **Figure 5B**). We evaluated the spike trains for these 20 PAs across all the five stimuli (**Figure 4**) by means of Victor-Purpura distance (VPd) (Victor and Purpura, 1996) spike train metrics. VPd gives a measure of the similarity between two spike trains by measuring the minimum cost necessary to transform one spike pattern into the other. The distance calculation is based on two rules, (i) adding/deleting a spike (cost = 1) and; (ii) shifting the spike by an interval of 1t (cost = q. 1t, in our simulation q = 10/s). By making comparisons between all the primary afferent spike trains evoked by different stimuli, we obtained a matrix of comparisons between all of the stimuli used (Rongala et al., 2017). For each neuron configuration, these matrices were used to compute the correlation index value for sensors with HEW synapses (**Figure 5C**, green) as given by computing the average VPd across all the 10 synapses with HEW for all the 5 stimuli. The same procedure was carried out to calculate the correlation index value for sensors with LEW synapses (**Figure 5C**, red). The difference between correlation index values of LEW and HEW sensors was the correlation index difference (**Figure 13**).

# Correlations Between Synaptic Responses

The correlation measure (**Figures 7**, **8**) was used to compute how correlated the synaptic responses were across different learning conditions. To compute this measure, we consider all the 80 synapses. We convolute the spike trains of all 80 PAs using the same EPSP time constants as in CN in vivo recordings (Bengtsson et al., 2013). Further we multiply each of these convoluted signals with the respective synaptic weights of the corresponding PA. The sum of all the 80 weighted convoluted signals was the simulated intracellular signal displayed (**Figures 7**, **8**). The correlation between two simulated intracellular signals is computed with the inbuilt MATLAB R function xcorr (with zero lag).

# Multi-Dimensional Scaling

To further illustrate the specific distribution of synaptic weights, we used the Classical multi-dimensional scaling (MDS) method (cmdscale, an inbuilt MATLAB R function). The distribution of weights (wexc) across the 80 PA synapses is here denoted the synaptic weight landscape. Using MDS, the differences between the synaptic landscapes of two or more CNs were visualized as distances in two-dimensional displays (**Figures 9A,C**,**10A**). The input distance vector to multi-dimensional scaling is calculated as a simple Euclidian distance between synaptic weight landscapes.

### Computation of the Sparseness of the Synaptic Weight Distribution

We also measured the degree of dispersion of synaptic weights, or the sparseness of the synaptic weights, among the 80 available PA synapses (**Figures 9B,D**, **10B**). For this purpose, we took the measured synaptic weight sparseness using the ratio between the l2 and l1 norm of the weight vector (Yin et al., 2014). This measure will report its minimal value when all synapses have exactly the same weight and a value of one when one synapse has maximum weight and all other synapses have zero weight. It was used to track the evolution of the synaptic weight dispersion during consecutive stimulus presentations and was therefore also an indicator of the stability of the learning.

#### Statistical Tests

After the learning process, the correlation index of the High End Weight (HEW) PA synapses (N = 10) was compared to that of the Low End Weight (LEW) PA synapses (N = 10) for each CN configuration using paired t-test (these distributions were on visual inspection comparable to normal distributions, see **Figure 5D**). The changes in correlation between the synaptic responses, of different CNs to the same stimulus and of the same CN to different stimuli, induced by training were statistically quantified using the paired, two-sided Wilcoxon signed rank test.

# RESULTS

### Inhibitory Synaptic Inputs and Intrinsic Responses of Cuneate Neurons in Vivo

To support the construction of our model, we started out by extending a previous characterization of the projection neurons of the cuneate nucleus in the decerebrate cat (Bengtsson et al., 2013) to also include their inhibitory synaptic inputs and intrinsic membrane responses, using a set of in vivo whole cell intracellular recordings from projection neurons (N = 15) and interneurons (N = 8) (**Figure 1A**). A touch to the receptive field of a cuneate projection neuron activated both excitatory and inhibitory synaptic inputs (**Figures 1A,B**) and the interneurons, which are responsible for synaptic inhibition evoked by skin stimulation in this preparation (Bengtsson et al., 2013). Spontaneous unitary inhibitory postsynaptic potentials (IPSPs), presumably driven by spontaneous interneuron spiking, were very small (**Figures 1B,C**) (the average median IPSP amplitude was −0.55 ± 0.08 mV, recorded from N = 15 cuneate projection neurons). In fact, this amplitude is an overestimate since many of the apparent IPSPs were too small to be detected from baseline noise using our template based identification method. In contrast, maximal evoked compound IPSPs (evoked from outside the receptive field of the excitatory input, see Bengtsson et al., 2013) were large (−14.3 ± 2.1 mV) and recruited in a continuous, gradable fashion (**Figure 1D**). Such maximal compound IPSPs had a response onset latency of 6.1 ± 0.7 ms, which is about 1.0 ms longer than the latency time of primary afferent EPSPs in the same type of preparation (Bengtsson et al., 2013). This observation strongly suggests that the compound IPSPs were due to primary afferent activation of the inhibitory interneurons. Together, these findings indicated that a high number of interneurons (>30, but possibly 100's) innervated each projection neuron with relatively uniformly low weighted inhibitory synapses.

A strong involvement of intrinsic voltage-activated conductances in shaping the membrane potential trajectory over time was suggested by the subthreshold spontaneous activity of the cuneate projection neurons recorded near resting potential (**Figure 1E**). Note the doublets, triplets or quadruplets of spike firing, which were observed in all projection neurons, but none of the interneurons. Occasional single spike firing events dramatically altered the spike afterhyperpolarization phase (red trace in **Figure 1E**), which further suggested the presence of intrinsic membrane potential dynamics. Indeed, using current clamp to generate temporary hyperpolarization of the membrane potential resulted in strong rebound excitation (**Figure 1E**), which demonstrates the contribution of such active intrinsic conductances. Similar intrinsic responses were recorded in all cuneate projection neurons when hyperpolarized by at least −200 pA for at least 200 ms (N = 15). Since such intrinsic responses can be expected to be an important part of defining the time course of calcium influx into the neuron, which in turn will be important for its learning process during a time-varying synaptic input, an important component of our full model below was the "intermediate dynamics" model of the of the cuneate projection neurons. We observed that there were relatively wide differences in the intensity and the speed of the intrinsic responses (**Figure 1E**) between the cuneate projection neurons. Therefore, we designed the intermediate dynamics model of these membrane potential dynamics against one of our recorded cuneate projection neurons, with the primary aim of approximating its overall dynamics across a range of conditions rather than focusing on capturing its behavior under any single condition exactly (for the full model, we later also will show model experiments in which parameters influencing the intermediate dynamics are adjusted to modify the intrinsic responses of the cuneate projection neurons, showing that model behavior was robust across a range of intermediate dynamics). Our intermediate dynamics model approximated the observed membrane potential dynamics well across a range of hyperpolarizations (**Figure 1F**). The model also resulted in a reproduction at least of the overall dynamics for responses evoked by touch (Figure S1).

In order to verify that the recorded IPSPs (**Figure 1**) could be generated by the local interneurons, we made intracellular recordings also from these neurons. Intracellular recordings were identified as projection neurons and interneurons based on their characteristic spike shapes (**Figure 2A**). The spontaneous spike firing frequencies were 11 ± 4.8 Hz (mean ± sd.) vs. 8.8 ± 2.2 Hz for projection neurons (N = 15) and interneurons (N = 8), respectively. Interneurons responded well to touch to the skin (**Figure 2B**) and on electrical skin stimulation (**Figure 2C**). The EPSPs evoked by electrical skin stimulation in the interneurons were considered monosynaptic based on the response latency time and the low variability of the EPSPs (Bengtsson et al., 2013; **Figure 2C**). Based on a qualitative analysis of the extent of their receptive fields on the skin (**Figure 2B**), interneurons were found to have much wider excitatory receptive fields than the cuneate projection neurons (which were previously quantitatively defined in Bengtsson et al., 2013). Combined with the rapid responses of the interneurons to PA synaptic input (**Figure 2C**), these observations made it plausible that their activation was responsible for the compound IPSPs recorded in the cuneate projection neurons on stimulation within or in the vicinity of their receptive fields (**Figure 1D**). The apparent absence of intrinsic regenerative responses to skin stimulation

FIGURE 1 | Inhibitory synaptic inputs and intrinsic responses of cuneate neurons in vivo. (A) Schematic of cuneate nucleus neuronal circuitry. Lines from the left symbolize axons originating from skin sensors. Their synaptic weights are indicated by the size of the triangles or as a stop ending for near-zero weights. Inhibitory interneurons (red for highlighted neuron, the others are gray) send axons that make inhibitory synapses on the projection neuron (black triangles). In the middle is a cuneate projection neuron, which sends its axon to the thalamus or the cerebellum. Inset: 3D illustration of the lower brainstem and the location of the cuneate nucleus (green). 3D scale bar indicate distances (1 mm per division). Yellow is gracile nucleus, small red volume is external cuneate nucleus. (B) Response of a (Continued) FIGURE 1 | projection neuron in vivo to a light skin touch within its receptive area. Arrows indicate putative spontaneous IPSPs, which are so small they are barely visible at this magnification of the voltage trace. (C) At a different magnification, examples of spontaneous IPSPs from one cuneate projection neuron are superimposed and averaged to the left. The peak amplitudes of 500 consecutive spontaneous IPSPs are shown in the box plot at right. (D) The gradual recruitment of summated, or compound, IPSPs with increased electrical stimulation intensity to a skin area adjacent to the excitatory receptive skin area of a sample projection neuron (black curve). Gray curve illustrates corresponding recruitment of unitary EPSPs (Bengtsson et al., 2013). Inset traces, average evoked IPSPs (averages of 20-50 repetitions) at different stimulation intensities. (E) Spontaneous activity of a projection neuron with zoomed-in, superimposed spikes to the right. Red trace indicates a case of a spontaneous single spike event, whereas most spontaneous firing occurred in doublets with an associated difference in voltage trajectory after the spike. Bottom, examples of rebound responses elicited by release from hyperpolarizing current injections (−400 pA for 200 ms). (F) Comparison between the rebound responses of a cuneate neuron in vivo and the CN model.

in the interneurons (**Figure 2B**) suggested the function of the cuneate interneurons to be a linear transfer of PA input, inverted to IPSPs, to the projection neurons. These findings motivated our simplification of the inhibitory synaptic inputs to the cuneate projection neurons to be modeled as a set of unitary weight inhibitory synaptic inputs directly generated by the spikes of the individual PAs (**Figures 3A,B**).

# Functional Structure of the CN Model

In our integrated cuneate model, the biological sensorized skin was replaced with a bionic fingertip (Oddo et al., 2016; **Figure 3A** and Figure S2). The model cuneate projection neuron (CN) was composed of the intermediate dynamics model described above (**Figure 1F**), that implemented input- and time-varying intrinsic responsiveness, and the excitatory and inhibitory synaptic inputs from the sensors that drove the activation of the intermediate dynamics of the CN ("Total calcium activity" in **Figure 3B**). The factor gmax allowed an adjustment of the ratio between the synaptic conductances and the background, leak, conductance, which is later in the paper used for simulating CNs with different input resistances or sizes. The rest of the diagram in **Figure 3B** indicates adaptation mechanisms for the excitatory and inhibitory synaptic weights, which implemented a Hebbiantype of learning combined with synaptic weight scaling as described below.

Depending on the correlation between the calcium activity in the main compartment and in the individual synaptic spaces (**Figures 3C,D**), the weight of each excitatory synapse was updated after each stimulus presentation (**Figures 3E–H**). To illustrate this process, **Figure 3E** shows that the local synaptic calcium activity crossed the learning threshold twice in this example. But for calculating the learning signal, also the calcium activity in the main compartment (**Figure 3D**) was taken into account. For the first episode of suprathreshold calcium activity in the synapse (**Figure 3E**), the calcium activity in the main compartment was in the "positive zone," i.e., sufficiently high to permit potentiation of the synaptic weight. Therefore, this first episode made the net learning drive of this synapse positive (**Figure 3F**). The second episode of suprathreshold calcium elevation in the synapse (**Figure 3E**) instead coincided with calcium activity in the main compartment being partly in the "negative zone" (**Figure 3D**). Therefore, the net learning drive for this episode was partly negative (**Figure 3F**). The net learning effect of each stimulus presentation and each synapse was the integral of all of these episodes, which in this case ended up positive (**Figure 3F**). For another synapse, the suprathreshold calcium activity (**Figure 3G**) mostly coincided with main compartment calcium activity being in the negative zone (**Figure 3D**), which resulted in a negative integral of the net learning drive (**Figure 3H**).

For each update of synaptic weights, the weight of the synapse itself was an additional factor used to scale the synaptic weight change (**Figure 3I**; synaptic weight compensation in **Figure 3B**) (Bi and Poo, 1998). Altogether, this architecture resulted in a Hebbian form of plasticity (Hebb, 1949). The end result was a form of plasticity that had properties being related to spike-time dependent plasticity (STDP) (Markram

FIGURE 3 | Functional structure of the CN model. (A) Structure of the CN network model. Colored lines indicate PA connections between the physical sensors of the bionic fingertip and the neurons of the model, blue triangles indicate variable weight synapses. Only cuneate projection neurons (CN1–CNn) were simulated individually, (Continued)

FIGURE 3 | where the number identifies a specific initial synaptic weight configuration (see Figure 5B). PA inputs were also provided to the CNs as a lump inhibitory synaptic input via an interneuron (red). (B) Graphic representation of the functional structure of the CN model. Components indicated with a dashed outline were the only elements with fixed parameters, whereas parameters with solid outlines were adjusted in the later part of this paper to simulate CNs with different intrinsic properties. (C) Subcellular factors involved in the synaptic plasticity of the CN, with variable weight PA synapses (blue) and inhibitory synapses (red). The CN neuron is divided into a main compartment, with reactive LVA and CAP conductances and synaptic spaces containing VGCCs (Higley and Sabatini, 2012) and a variable number of AMPAr:s. AMPAr, excitatory glutamate receptors; VGCCs, voltage gated calcium channels; GABAa, inhibitory synaptic receptors; LVA, low–threshold voltage activated calcium channels; CAP, calcium-dependent potassium channels. (D) The calcium activity of the main compartment varied over time due to synaptic input and the responses they elicited via the intermediate dynamics model of the CN. The level of the calcium activity, in relationship to the learning threshold, defined when the cell was in the positive zone (i.e., potentiation mode) or in the negative zone (depression mode). (E) The local calcium activity threshold used to define the eligibility for plasticity for each PA synapse. (F) The net learning drive for the PA synapse varied depending on the temporal correlation between the zero-offset main compartment calcium activity and the local, or synapse space, calcium activity. (G) The activity of another PA synapse under the same stimulus presentation and; (H) its net learning drive. (I) The synaptic weight compensation constant that was multiplied with the integral net learning drive to calculate the final weight change per stimulus presentation. For (D,E,G), note that the y axes indicate the relative magnitude of calcium signals for each compartment and not an actual estimation of the calcium concentration.

et al., 1997; Graupner and Brunel, 2012) but being less sharply dependent on the exact spike timing (Helias et al., 2008). The motivation for taking this approach to formulate the driving forces for the learning rule was to let the dynamics of the cellular learning signal, and its shaping by the integration of peripheral synaptic input have its origin in biological data recorded from cuneate projection neurons in vivo (**Figures 1E,F**) rather than in assumptions made for other learning systems. We believe that overall, net learning effects are comparable to those obtained under some settings of those previous models (Graupner and Brunel, 2012) (Discussion) but more specifically constrained to the effects that might arise in the cuneate projection neurons.

The integrated model also used two separate feedback regulatory systems. First, the total inhibitory synaptic weight was updated based on the total calcium activity, according to a fixed set point (**Figure 3B**). Adaptation of the total inhibitory synaptic weights was implemented because we found that the change in the weights of the excitatory PA synapses on the CN neuron could lead to large changes in the levels of activity in the voltage-activated calcium channels, which in turn led to unstable learning. A plastic regulation of the weight of inhibitory synapses against postsynaptic calcium activity or postsynaptic calcium-dependent bursting activity has previously been found in various neurons across the brain (Kawaguchi et al., 2011; Hulme and Connelly, 2014; Lourenço et al., 2014).

Secondly, SetPointP wexc is a control signal used to implement a homeostatic synaptic plasticity mechanisms (Turrigiano and Nelson, 2004; Turrigiano, 2011) through which the sum of the excitatory synaptic weights is kept under control. The total excitatory synaptic weight was driven toward this set point by adjusting the learning threshold for synaptic potentiation vs. depression (**Figure 3B**; positive vs. negative zone in **Figure 3D**). Hence, the propensity for LTP or LTD was affected by the current total synaptic weight and the effect of the learning threshold was therefore a synaptic weight scaling. The set point of the learning threshold (**Figure 3B**) is the fixed sum of the excitatory synaptic weights toward which the adaptation of the learning threshold strives to bring the system. Further details and motivations are given in Methods.

# Synaptic Weight Changes During the Learning Process

The sensorized bionic fingertip was used to provide the CN with training data. The bionic fingertip provided the important feature of a system of PA sensors where there is a consistent relationship between the activation of the different sensors across different conditions or tactile experiences. The CN learning process consisted of 1,500 presentations of five training stimuli, each activating all of the 80 PA sensory inputs but in different spatiotemporal patterns (**Figure 4**), presented in a pseudorandom order. The training stimuli consisted in three touch-an-slide conditions and in two dynamic indentation stimuli (**Figure 4**). During the learning process, the gradual transition from initial, or "seed," synaptic weights to the final distribution of synaptic weights ("end weight") was monitored (**Figure 5A**). Interestingly, the weight of a synapse could for example evolve from a very low seed weight to a high end weight (green thick line), or from high seed weight to nearzero end weight (red thick line) and any intermediate or mixed development. Such examples showed that the fate of an individual synaptic weight during the learning process was not depending on its starting point but suggested some underlying functional principle.

The learning process resulted in a transformation of the initial random distribution of low-moderate weight synapses (seed weight, SW) (**Figure 5B**) to a highly skewed distribution with a few high weight synapses (end weight, EW) (**Figure 5C**), similar to cuneate neurons in adult animals (Bengtsson et al., 2013). We simulated CNs with five different pseudo-randomized SWs, which below are referred to as CN1-5, based on the assumption that initial random synaptic weight distributions likely differ between biological cuneate neurons in the same animal. For each CN, there was at the end of the learning process a unique end weight distribution (EW), dominated by different specific synapses with high weights (Figure S3). The seed weight losers (red bars in **Figure 5B**) differed from the end weight winners (green bars in **Figure 5C** and Figure S3) in that they conveyed sensor spiking activity that was less correlated across the training stimuli (**Figure 5D**), as measured by the Victor-Purpura distance metric for comparisons between spike trains. During the learning process, the overall adaptive regulation of the total excitatory and inhibitory synaptic weights in the model

(**Figure 3B**) helped stabilizing the activity of the CN (**Figure 6**), which was important for a stable learning outcome. In addition, in parallel with the acquisition of highly skewed synaptic weight distributions, the cuneate neurons also gradually acquired a firing that was dominated by brief, episodic bursts of spike output, similar to the firing mode observed for these neurons in adulthood in vivo (Bengtsson et al., 2013; **Figures 1**, **6D**, Figure S1).

## Learning Resulted in Decorrelations of Synaptic Inputs Between Neurons and Stimuli

The next question asked was what functional effects the learning process resulted in. We found that between the individual CNs, the learning process resulted in a decorrelation of the initially highly correlated temporal patterns of the synaptic responses to each given training stimulus (**Figures 7A,B**).

before and after learning in the example CN (CN1).

Hence, a consequence of the learning process was an increased differentiation of the CN responses to a given input, which was statistically significant at p ≤ 0.001 (Wilcoxon signed rank test). In addition to testing the stimuli that the CNs were trained to, we also used a set of previously non-encountered stimuli, consisting of three types of texture sliding and two touch indentations, here referred to as non-training stimuli. Similar to the case with the training stimuli, the learning process resulted in a decorrelation of the synaptic responses of the CNs also to the non-training stimuli (**Figures 7C,D**; p < 0.001, Wilcoxon signed rank test). Moreover, in the individual CN, whereas the differences in the spatiotemporal patterns of PA input between stimuli to some extent generated specific synaptic responses already in the "naïve" state, the responses became markedly differentiated with learning (**Figures 8A,B**). For each of the five CNs, CN1-5, the learning process resulted in a statistically significant decrease in the correlation between the synaptic responses to the five training stimuli (p < 0.001, Wilcoxon signed rank test). Such synaptic response differentiation is a basis for tactile input feature segregation (Jörntell et al., 2014)