Learning to Estimate Dynamical State with Probabilistic Population Codes
Identifieur interne : 003D39 ( Ncbi/Merge ); précédent : 003D38; suivant : 003D40Learning to Estimate Dynamical State with Probabilistic Population Codes
Auteurs : Joseph G. Makin [États-Unis] ; Benjamin K. Dichter [États-Unis] ; Philip N. Sabes [États-Unis]Source :
- PLoS Computational Biology [ 1553-734X ] ; 2015.
Abstract
Tracking moving objects, including one’s own body, is a fundamental ability of higher organisms, playing a central role in many perceptual and motor tasks. While it is unknown how the brain learns to follow and predict the dynamics of objects, it is known that this process of state estimation can be learned purely from the statistics of noisy observations. When the dynamics are simply linear with additive Gaussian noise, the optimal solution is the well known Kalman filter (KF), the parameters of which can be learned via latent-variable density estimation (the EM algorithm). The brain does not, however, directly manipulate matrices and vectors, but instead appears to represent probability distributions with the firing rates of population of neurons, “probabilistic population codes.” We show that a recurrent neural network—a modified form of an exponential family harmonium (EFH)—that takes a linear probabilistic population code as input can learn, without supervision, to estimate the state of a linear dynamical system. After observing a series of population responses (spike counts) to the position of a moving object, the network learns to represent the velocity of the object and forms nearly optimal predictions about the position at the next time-step. This result builds on our previous work showing that a similar network can learn to perform multisensory integration and coordinate transformations for static stimuli. The receptive fields of the trained network also make qualitative predictions about the developing and learning brain: tuning gradually emerges for higher-order dynamical states not explicitly present in the inputs, appearing as
Url:
DOI: 10.1371/journal.pcbi.1004554
PubMed: 26540152
PubMed Central: 4634970
Links toward previous steps (curation, corpus...)
- to stream Pmc, to step Corpus: 000262
- to stream Pmc, to step Curation: 000262
- to stream Pmc, to step Checkpoint: 000476
Links to Exploration step
PMC:4634970Le document en format XML
<record><TEI><teiHeader><fileDesc><titleStmt><title xml:lang="en">Learning to Estimate Dynamical State with Probabilistic Population Codes</title>
<author><name sortKey="Makin, Joseph G" sort="Makin, Joseph G" uniqKey="Makin J" first="Joseph G." last="Makin">Joseph G. Makin</name>
<affiliation wicri:level="2"><nlm:aff id="aff001"><addr-line>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff002"><addr-line>Department of Physiology, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Department of Physiology, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
</author>
<author><name sortKey="Dichter, Benjamin K" sort="Dichter, Benjamin K" uniqKey="Dichter B" first="Benjamin K." last="Dichter">Benjamin K. Dichter</name>
<affiliation wicri:level="2"><nlm:aff id="aff001"><addr-line>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff003"><addr-line>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
</author>
<author><name sortKey="Sabes, Philip N" sort="Sabes, Philip N" uniqKey="Sabes P" first="Philip N." last="Sabes">Philip N. Sabes</name>
<affiliation wicri:level="2"><nlm:aff id="aff001"><addr-line>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff002"><addr-line>Department of Physiology, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Department of Physiology, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff003"><addr-line>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
</author>
</titleStmt>
<publicationStmt><idno type="wicri:source">PMC</idno>
<idno type="pmid">26540152</idno>
<idno type="pmc">4634970</idno>
<idno type="url">http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4634970</idno>
<idno type="RBID">PMC:4634970</idno>
<idno type="doi">10.1371/journal.pcbi.1004554</idno>
<date when="2015">2015</date>
<idno type="wicri:Area/Pmc/Corpus">000262</idno>
<idno type="wicri:Area/Pmc/Curation">000262</idno>
<idno type="wicri:Area/Pmc/Checkpoint">000476</idno>
<idno type="wicri:Area/Ncbi/Merge">003D39</idno>
</publicationStmt>
<sourceDesc><biblStruct><analytic><title xml:lang="en" level="a" type="main">Learning to Estimate Dynamical State with Probabilistic Population Codes</title>
<author><name sortKey="Makin, Joseph G" sort="Makin, Joseph G" uniqKey="Makin J" first="Joseph G." last="Makin">Joseph G. Makin</name>
<affiliation wicri:level="2"><nlm:aff id="aff001"><addr-line>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff002"><addr-line>Department of Physiology, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Department of Physiology, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
</author>
<author><name sortKey="Dichter, Benjamin K" sort="Dichter, Benjamin K" uniqKey="Dichter B" first="Benjamin K." last="Dichter">Benjamin K. Dichter</name>
<affiliation wicri:level="2"><nlm:aff id="aff001"><addr-line>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff003"><addr-line>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
</author>
<author><name sortKey="Sabes, Philip N" sort="Sabes, Philip N" uniqKey="Sabes P" first="Philip N." last="Sabes">Philip N. Sabes</name>
<affiliation wicri:level="2"><nlm:aff id="aff001"><addr-line>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff002"><addr-line>Department of Physiology, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>Department of Physiology, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
<affiliation wicri:level="2"><nlm:aff id="aff003"><addr-line>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</nlm:aff>
<country xml:lang="fr">États-Unis</country>
<wicri:regionArea>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California</wicri:regionArea>
<placeName><region type="state">Californie</region>
</placeName>
</affiliation>
</author>
</analytic>
<series><title level="j">PLoS Computational Biology</title>
<idno type="ISSN">1553-734X</idno>
<idno type="eISSN">1553-7358</idno>
<imprint><date when="2015">2015</date>
</imprint>
</series>
</biblStruct>
</sourceDesc>
</fileDesc>
<profileDesc><textClass></textClass>
</profileDesc>
</teiHeader>
<front><div type="abstract" xml:lang="en"><p>Tracking moving objects, including one’s own body, is a fundamental ability of higher organisms, playing a central role in many perceptual and motor tasks. While it is unknown how the brain learns to follow and predict the dynamics of objects, it is known that this process of state estimation can be learned purely from the statistics of noisy observations. When the dynamics are simply linear with additive Gaussian noise, the optimal solution is the well known Kalman filter (KF), the parameters of which can be learned via latent-variable density estimation (the EM algorithm). The brain does not, however, directly manipulate matrices and vectors, but instead appears to represent probability distributions with the firing rates of population of neurons, “probabilistic population codes.” We show that a recurrent neural network—a modified form of an exponential family harmonium (EFH)—that takes a linear probabilistic population code as input can learn, without supervision, to estimate the state of a linear dynamical system. After observing a series of population responses (spike counts) to the position of a moving object, the network learns to represent the velocity of the object and forms nearly optimal predictions about the position at the next time-step. This result builds on our previous work showing that a similar network can learn to perform multisensory integration and coordinate transformations for static stimuli. The receptive fields of the trained network also make qualitative predictions about the developing and learning brain: tuning gradually emerges for higher-order dynamical states not explicitly present in the inputs, appearing as <italic>delayed</italic>
tuning for the lower-order states.</p>
</div>
</front>
<back><div1 type="bibliography"><listBibl><biblStruct><analytic><author><name sortKey="Foldiak, P" uniqKey="Foldiak P">P Földiák</name>
</author>
<author><name sortKey="Eeckman, Fh" uniqKey="Eeckman F">FH Eeckman</name>
</author>
<author><name sortKey="Bower, Jm" uniqKey="Bower J">JM Bower</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Ma, Wj" uniqKey="Ma W">WJ Ma</name>
</author>
<author><name sortKey="Beck, Jm" uniqKey="Beck J">JM Beck</name>
</author>
<author><name sortKey="Latham, Pe" uniqKey="Latham P">PE Latham</name>
</author>
<author><name sortKey="Pouget, A" uniqKey="Pouget A">A Pouget</name>
</author>
</analytic>
</biblStruct>
<biblStruct></biblStruct>
<biblStruct><analytic><author><name sortKey="Makin, Jg" uniqKey="Makin J">JG Makin</name>
</author>
<author><name sortKey="Fellows, Mr" uniqKey="Fellows M">MR Fellows</name>
</author>
<author><name sortKey="Sabes, Pn" uniqKey="Sabes P">PN Sabes</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Van Beers, Rj" uniqKey="Van Beers R">RJ van Beers</name>
</author>
<author><name sortKey="Sittig, A" uniqKey="Sittig A">A Sittig</name>
</author>
<author><name sortKey="Van Der Gon, Jjd" uniqKey="Van Der Gon J">JJD van Der Gon</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Ernst, Mo" uniqKey="Ernst M">MO Ernst</name>
</author>
<author><name sortKey="Banks, Ms" uniqKey="Banks M">MS Banks</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Beck, Jm" uniqKey="Beck J">JM Beck</name>
</author>
<author><name sortKey="Latham, Pe" uniqKey="Latham P">PE Latham</name>
</author>
<author><name sortKey="Pouget, A" uniqKey="Pouget A">A Pouget</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Boerlin, M" uniqKey="Boerlin M">M Boerlin</name>
</author>
<author><name sortKey="Deneve, S" uniqKey="Deneve S">S Denève</name>
</author>
</analytic>
</biblStruct>
<biblStruct></biblStruct>
<biblStruct><analytic><author><name sortKey="Wolpert, Dm" uniqKey="Wolpert D">DM Wolpert</name>
</author>
<author><name sortKey="Goodbody, Sj" uniqKey="Goodbody S">SJ Goodbody</name>
</author>
<author><name sortKey="Husain, M" uniqKey="Husain M">M Husain</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Mulliken, Gh" uniqKey="Mulliken G">GH Mulliken</name>
</author>
<author><name sortKey="Musallam, S" uniqKey="Musallam S">S Musallam</name>
</author>
<author><name sortKey="Andersen, Ra" uniqKey="Andersen R">RA Andersen</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Hinton, Ge" uniqKey="Hinton G">GE Hinton</name>
</author>
<author><name sortKey="Osindero, S" uniqKey="Osindero S">S Osindero</name>
</author>
<author><name sortKey="Teh, Yw" uniqKey="Teh Y">YW Teh</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Hinton, Ge" uniqKey="Hinton G">GE Hinton</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Scherpen, Jma" uniqKey="Scherpen J">JMA Scherpen</name>
</author>
<author><name sortKey="Levine, Ws" uniqKey="Levine W">WS Levine</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Hinton, Ge" uniqKey="Hinton G">GE Hinton</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Makin, Jg" uniqKey="Makin J">JG Makin</name>
</author>
<author><name sortKey="Fellows, Mr" uniqKey="Fellows M">MR Fellows</name>
</author>
<author><name sortKey="Sabes, Pn" uniqKey="Sabes P">PN Sabes</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Dayan, P" uniqKey="Dayan P">P Dayan</name>
</author>
<author><name sortKey="Abbott, Lf" uniqKey="Abbott L">LF Abbott</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Ghahramani, Z" uniqKey="Ghahramani Z">Z Ghahramani</name>
</author>
<author><name sortKey="Hinton, Ge" uniqKey="Hinton G">GE Hinton</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Cover, Tm" uniqKey="Cover T">TM Cover</name>
</author>
<author><name sortKey="Thomas, Ja" uniqKey="Thomas J">JA Thomas</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Zar, Jh" uniqKey="Zar J">JH Zar</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Andersen, Ra" uniqKey="Andersen R">RA Andersen</name>
</author>
<author><name sortKey="Snyder, Lh" uniqKey="Snyder L">LH Snyder</name>
</author>
<author><name sortKey="Bradley, Dc" uniqKey="Bradley D">DC Bradley</name>
</author>
<author><name sortKey="Xing, J" uniqKey="Xing J">J Xing</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Egger, Sw" uniqKey="Egger S">SW Egger</name>
</author>
<author><name sortKey="Britten, Kh" uniqKey="Britten K">KH Britten</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Kuenzle, H" uniqKey="Kuenzle H">H Kuenzle</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Ghosh, S" uniqKey="Ghosh S">S Ghosh</name>
</author>
<author><name sortKey="Brinkman, C" uniqKey="Brinkman C">C Brinkman</name>
</author>
<author><name sortKey="Porter, R" uniqKey="Porter R">R Porter</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Burchinskaya, Lf" uniqKey="Burchinskaya L">LF Burchinskaya</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Carracedo, Lm" uniqKey="Carracedo L">LM Carracedo</name>
</author>
<author><name sortKey="Kjeldsen, H" uniqKey="Kjeldsen H">H Kjeldsen</name>
</author>
<author><name sortKey="Cunnington, L" uniqKey="Cunnington L">L Cunnington</name>
</author>
<author><name sortKey="Jenkins, A" uniqKey="Jenkins A">a Jenkins</name>
</author>
<author><name sortKey="Schofield, I" uniqKey="Schofield I">I Schofield</name>
</author>
<author><name sortKey="Cunningham, Mo" uniqKey="Cunningham M">MO Cunningham</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Rao, Rpn" uniqKey="Rao R">RPN Rao</name>
</author>
<author><name sortKey="Ballard, Dh" uniqKey="Ballard D">DH Ballard</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Deneve, S" uniqKey="Deneve S">S Denève</name>
</author>
<author><name sortKey="Duhamel, Jr" uniqKey="Duhamel J">JR Duhamel</name>
</author>
<author><name sortKey="Pouget, A" uniqKey="Pouget A">A Pouget</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Huys, Qjm" uniqKey="Huys Q">QJM Huys</name>
</author>
<author><name sortKey="Zemel, Rs" uniqKey="Zemel R">RS Zemel</name>
</author>
<author><name sortKey="Natarajan, R" uniqKey="Natarajan R">R Natarajan</name>
</author>
<author><name sortKey="Dayan, P" uniqKey="Dayan P">P Dayan</name>
</author>
</analytic>
</biblStruct>
<biblStruct><analytic><author><name sortKey="Natarajan, R" uniqKey="Natarajan R">R Natarajan</name>
</author>
<author><name sortKey="Huys, Qjm" uniqKey="Huys Q">QJM Huys</name>
</author>
<author><name sortKey="Dayan, P" uniqKey="Dayan P">P Dayan</name>
</author>
<author><name sortKey="Zemel, Rs" uniqKey="Zemel R">RS Zemel</name>
</author>
</analytic>
</biblStruct>
<biblStruct></biblStruct>
<biblStruct><analytic><author><name sortKey="Sutskever, I" uniqKey="Sutskever I">I Sutskever</name>
</author>
<author><name sortKey="Hinton, Ge" uniqKey="Hinton G">GE Hinton</name>
</author>
</analytic>
</biblStruct>
<biblStruct></biblStruct>
<biblStruct><analytic><author><name sortKey="Wolpert, Dm" uniqKey="Wolpert D">DM Wolpert</name>
</author>
<author><name sortKey="Ghahramani, Z" uniqKey="Ghahramani Z">Z Ghahramani</name>
</author>
<author><name sortKey="Jordan, Mi" uniqKey="Jordan M">MI Jordan</name>
</author>
</analytic>
</biblStruct>
</listBibl>
</div1>
</back>
</TEI>
<pmc article-type="research-article"><pmc-dir>properties open_access</pmc-dir>
<front><journal-meta><journal-id journal-id-type="nlm-ta">PLoS Comput Biol</journal-id>
<journal-id journal-id-type="iso-abbrev">PLoS Comput. Biol</journal-id>
<journal-id journal-id-type="publisher-id">plos</journal-id>
<journal-id journal-id-type="pmc">ploscomp</journal-id>
<journal-title-group><journal-title>PLoS Computational Biology</journal-title>
</journal-title-group>
<issn pub-type="ppub">1553-734X</issn>
<issn pub-type="epub">1553-7358</issn>
<publisher><publisher-name>Public Library of Science</publisher-name>
<publisher-loc>San Francisco, CA USA</publisher-loc>
</publisher>
</journal-meta>
<article-meta><article-id pub-id-type="pmid">26540152</article-id>
<article-id pub-id-type="pmc">4634970</article-id>
<article-id pub-id-type="publisher-id">PCOMPBIOL-D-15-00643</article-id>
<article-id pub-id-type="doi">10.1371/journal.pcbi.1004554</article-id>
<article-categories><subj-group subj-group-type="heading"><subject>Research Article</subject>
</subj-group>
</article-categories>
<title-group><article-title>Learning to Estimate Dynamical State with Probabilistic Population Codes</article-title>
<alt-title alt-title-type="running-head">Learning to Track Moving Stimuli with Population Codes</alt-title>
</title-group>
<contrib-group><contrib contrib-type="author" equal-contrib="yes"><name><surname>Makin</surname>
<given-names>Joseph G.</given-names>
</name>
<xref ref-type="corresp" rid="cor001">*</xref>
<xref ref-type="aff" rid="aff001"><sup>1</sup>
</xref>
<xref ref-type="aff" rid="aff002"><sup>2</sup>
</xref>
</contrib>
<contrib contrib-type="author" equal-contrib="yes"><name><surname>Dichter</surname>
<given-names>Benjamin K.</given-names>
</name>
<xref ref-type="aff" rid="aff001"><sup>1</sup>
</xref>
<xref ref-type="aff" rid="aff003"><sup>3</sup>
</xref>
</contrib>
<contrib contrib-type="author"><name><surname>Sabes</surname>
<given-names>Philip N.</given-names>
</name>
<xref ref-type="aff" rid="aff001"><sup>1</sup>
</xref>
<xref ref-type="aff" rid="aff002"><sup>2</sup>
</xref>
<xref ref-type="aff" rid="aff003"><sup>3</sup>
</xref>
</contrib>
</contrib-group>
<aff id="aff001"><label>1</label>
<addr-line>Center for Integrative Neuroscience, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</aff>
<aff id="aff002"><label>2</label>
<addr-line>Department of Physiology, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</aff>
<aff id="aff003"><label>3</label>
<addr-line>UC Berkeley-UCSF Graduate Program in Bioengineering, University of California, San Francisco, San Francisco, California, United States of America</addr-line>
</aff>
<contrib-group><contrib contrib-type="editor"><name><surname>Beck</surname>
<given-names>Jeff</given-names>
</name>
<role>Editor</role>
<xref ref-type="aff" rid="edit1"></xref>
</contrib>
</contrib-group>
<aff id="edit1"><addr-line>Duke University, UNITED STATES</addr-line>
</aff>
<author-notes><fn fn-type="conflict" id="coi001"><p>The authors have declared that no competing interests exist.</p>
</fn>
<fn fn-type="con" id="contrib001"><p>Conceived and designed the experiments: JGM BKD PNS. Performed the experiments: JGM BKD. Analyzed the data: JGM BKD PNS. Wrote the paper: JGM BKD PNS. Proposed the technique that extends the harmonium to dynamical data: BKD.</p>
</fn>
<corresp id="cor001">* E-mail: <email>makin@phy.ucsf.edu</email>
</corresp>
</author-notes>
<pub-date pub-type="collection"><month>11</month>
<year>2015</year>
</pub-date>
<pub-date pub-type="epub"><day>5</day>
<month>11</month>
<year>2015</year>
</pub-date>
<volume>11</volume>
<issue>11</issue>
<elocation-id>e1004554</elocation-id>
<history><date date-type="received"><day>20</day>
<month>4</month>
<year>2015</year>
</date>
<date date-type="accepted"><day>26</day>
<month>8</month>
<year>2015</year>
</date>
</history>
<permissions><copyright-year>2015</copyright-year>
<copyright-holder>Makin et al</copyright-holder>
<license xlink:href="http://creativecommons.org/licenses/by/4.0/"><license-p>This is an open access article distributed under the terms of the <ext-link ext-link-type="uri" xlink:href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution License</ext-link>
, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited</license-p>
</license>
</permissions>
<self-uri content-type="pdf" xlink:type="simple" xlink:href="pcbi.1004554.pdf"></self-uri>
<abstract><p>Tracking moving objects, including one’s own body, is a fundamental ability of higher organisms, playing a central role in many perceptual and motor tasks. While it is unknown how the brain learns to follow and predict the dynamics of objects, it is known that this process of state estimation can be learned purely from the statistics of noisy observations. When the dynamics are simply linear with additive Gaussian noise, the optimal solution is the well known Kalman filter (KF), the parameters of which can be learned via latent-variable density estimation (the EM algorithm). The brain does not, however, directly manipulate matrices and vectors, but instead appears to represent probability distributions with the firing rates of population of neurons, “probabilistic population codes.” We show that a recurrent neural network—a modified form of an exponential family harmonium (EFH)—that takes a linear probabilistic population code as input can learn, without supervision, to estimate the state of a linear dynamical system. After observing a series of population responses (spike counts) to the position of a moving object, the network learns to represent the velocity of the object and forms nearly optimal predictions about the position at the next time-step. This result builds on our previous work showing that a similar network can learn to perform multisensory integration and coordinate transformations for static stimuli. The receptive fields of the trained network also make qualitative predictions about the developing and learning brain: tuning gradually emerges for higher-order dynamical states not explicitly present in the inputs, appearing as <italic>delayed</italic>
tuning for the lower-order states.</p>
</abstract>
<abstract abstract-type="summary"><title>Author Summary</title>
<p>A basic task for animals is to track objects—predators, prey, even their own limbs—as they move through the world. Because the position estimates provided by the senses are not error-free, higher levels of performance can be, and are, achieved when the velocity and acceleration, as well as the position, of the object are taken into account. Likewise, tracking of limbs under voluntary control can be improved by considering the motor command that is (partially) responsible for its trajectory. Engineers have built tools to solve precisely these problems, and even to learn dynamical features of the object to be tracked. How does the brain do it? We show how artificial networks of neurons can learn to solve this task, simply by trying to become good predictive models of their incoming data—as long as some of those data are the activities of the neurons themselves at a fixed time delay, while the remainder (imperfectly) report the current position. The tracking scheme the network learns to use—keeping track of past positions; the corresponding receptive fields; and the manner in which they are learned, provide predictions for brain areas involved in tracking, like the posterior parietal cortex.</p>
</abstract>
<funding-group><funding-statement>This work was funded by: Defense Advanced Research Projects Agency (N66001-10-C-2010, W911NF-14-2-0043); National Institutes of Health Vision Training Grant (5T32EY007120-22); and UCSF Wheeler Center for the Neurobiology of Addiction (<ext-link ext-link-type="uri" xlink:href="http://physio.ucsf.edu/WheelerCenter/Center/">http://physio.ucsf.edu/WheelerCenter/Center/</ext-link>
). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.</funding-statement>
</funding-group>
<counts><fig-count count="9"></fig-count>
<table-count count="0"></table-count>
<page-count count="28"></page-count>
</counts>
<custom-meta-group><custom-meta id="data-availability"><meta-name>Data Availability</meta-name>
<meta-value>The code that generates all data, as well as trains and tests the models and produces the figures, is fully available at: <ext-link ext-link-type="uri" xlink:href="https://github.com/jgmakin/rbmish">https://github.com/jgmakin/rbmish</ext-link>
.</meta-value>
</custom-meta>
</custom-meta-group>
</article-meta>
<notes><title>Data Availability</title>
<p>The code that generates all data, as well as trains and tests the models and produces the figures, is fully available at: <ext-link ext-link-type="uri" xlink:href="https://github.com/jgmakin/rbmish">https://github.com/jgmakin/rbmish</ext-link>
.</p>
</notes>
</front>
<body><sec sec-type="intro" id="sec001"><title>Introduction</title>
<p>Over the last decade, neuroscience has come increasingly to believe that sensory systems represent not merely stimuli, but probability distributions over them. This conclusion follows from two observations. The first is that the apparent stochasticity of the response, <bold>R</bold>
, of a population of neurons inherently represents the likelihood of the stimulus <bold>s</bold>
: <bold>R</bold>
∼ <italic>p</italic>
(<bold>r</bold>
∣<bold>s</bold>
) [<xref rid="pcbi.1004554.ref001" ref-type="bibr">1</xref>
]. The second is that certain common computations essential to the function of many animals require keeping track of probability distributions over stimuli, rather than mere point estimates. For example, primates integrate information from multiple senses by weighting each sense by its reliability (inverse variance) [<xref rid="pcbi.1004554.ref005" ref-type="bibr">5</xref>
, <xref rid="pcbi.1004554.ref006" ref-type="bibr">6</xref>
]. This framework has been used to hand-wire neural networks that integrate spatial information across sensory modalities and across time [<xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
, <xref rid="pcbi.1004554.ref007" ref-type="bibr">7</xref>
, <xref rid="pcbi.1004554.ref008" ref-type="bibr">8</xref>
]. The more challenging problem faced by the brain, however, is to <italic>learn</italic>
to perform these tasks.</p>
<p>We have recently shown [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
, <xref rid="pcbi.1004554.ref009" ref-type="bibr">9</xref>
] that the problem of learning to integrate information about a common stimulus from multiple, unisensory populations of neurons can be solved by a neural network that implements a form of unsupervised learning called density estimation. Such a network learns to represent the joint probability density of the unisensory responses—to build a good model for these data—in terms of the activities of its downstream, multisensory units. For example [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
], an exponential family harmonium (EFH) [<xref rid="pcbi.1004554.ref003" ref-type="bibr">3</xref>
] trained on the activities of two populations of Gaussian-tuned, Poisson neurons (linear probabilistic population codes [<xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
]) that tile their respective sensory spaces (visual and proprioceptive, e.g.) will learn to extract the “common cause” of these populations, encoding the stimulus in its hidden layer. In this case, the unisensory information available on a “trial” can be characterized by two means (best estimates) and two variances (inverse reliabilities); and the estimate extracted by the hidden units of the trained network is precisely the inverse-variance-weighted convex combination that primates appear in psychophysical studies to use.</p>
<p>Ecologically, however, the critical challenge is not typically to estimate the location of a static object, but to track the state of a dynamically changing environment. This task likewise requires reliability-weighted combination of information, in this case of the current sensory evidence and the current best estimate of the state given past information. But it is considerably more difficult, since its solution requires learning a predictive model of the dynamics, which is not explicitly encoded in the sensory reports. In the case of Gaussian noise and linear dynamics (LDS), this recursive process is described by the Kalman filter, the parameters of which can be acquired with well-known iterative learning schemes. How the brain learns to solve this problem, however, is unknown.</p>
<p>Here we propose a neural model that accomplishes this task. We show that by adding recurrent connections to an EFH similar to that used in [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
], the network can learn to estimate the state of a dynamical system. For concreteness, we consider the problem of tracking the dynamical state of the upper limb, a necessary computation for accurate and precise movement planning and control. In this case, the neural circuit corresponds to the posterior parietal cortex (PPC), which appears to subserve state estimation [<xref rid="pcbi.1004554.ref010" ref-type="bibr">10</xref>
, <xref rid="pcbi.1004554.ref011" ref-type="bibr">11</xref>
]; and its inputs are taken to be a population of proprioceptive neurons. The network’s performance can be quantified precisely by restricting our view to linear-Gaussian dynamics, where the filtering and learning problems have known optimal solutions (respectively, the Kalman filter and expectation-maximization, a maximum-likelihood algorithm). And indeed, performance approaches that optimum.</p>
<p>We then extend the network to controlled dynamical systems. Under the assumption that the controls are provided by motor cortex, these too are observed only noisily by PPC, in the form of efference copy, which the network must then learn to interpret as motor commands. State estimation is again close to optimal. In addition, the network is neurally plausible in both its representation of stimulus probabilities [<xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
] and in the unsupervised learning procedure, which relies only on pairwise correlations between firing rates of connected neurons [<xref rid="pcbi.1004554.ref012" ref-type="bibr">12</xref>
, <xref rid="pcbi.1004554.ref013" ref-type="bibr">13</xref>
]. Finally, the network makes two predictions about neural circuits that learn to perform state estimation: (1) During learning, position receptive fields will emerge before velocity receptive fields; or more generally, receptive fields will develop from lower- to higher-order states, especially when explicit information about the higher-order states is not in the inputs. (2) Filtering is implemented by tuning to past positions (or more generally, lower-order states), rather than tuning directly to velocity (or more generally, higher-order states).</p>
</sec>
<sec sec-type="results" id="sec002"><title>Results</title>
<sec id="sec003"><title>Network performance</title>
<sec id="sec004"><title>The filtering problem</title>
<p>We present results for an uncontrolled and a controlled dynamical system. For both, the basic dynamical system is second-order (position, velocity), discrete-time, stochastic, and linear. The noise in state transitions is additive, white, and Gaussian. The “observation” at time <italic>t</italic>
is the response (<inline-formula id="pcbi.1004554.e001"><alternatives><graphic xlink:href="pcbi.1004554.e001.jpg" id="pcbi.1004554.e001g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M1"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) of a population of Poisson neurons, with Gaussian (bell-shaped) tuning curves that smoothly tile position. Since these neurons are taken to be reporting the proprioceptive sense, the “position” variable is the angle of a (single) joint, Θ. For the controlled system, there is a second population of Poisson neurons—carrying the “efference copy,” <inline-formula id="pcbi.1004554.e002"><alternatives><graphic xlink:href="pcbi.1004554.e002.jpg" id="pcbi.1004554.e002g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M2"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
—that smoothly tiles a space of input torques, <bold>U</bold>
. Details appear in the methods section <bold>Input-data generation</bold>
.</p>
<p>The task of tracking an object (<italic>estimation</italic>
) is to provide, at time <italic>t</italic>
, the best estimate of the location that can be computed from all the noisy observations from time 0 up to <italic>t</italic>
. When current position depends on only a finite number <italic>n</italic>
of past positions, this problem can be solved recursively: rather than retaining a history of all past observations, it is necessary to maintain only the current best estimate of the state (a vector whose dimension is set by <italic>n</italic>
), and the reliabilities of these estimates (a covariance matrix). By artful design of the system (see <xref ref-type="sec" rid="sec004">Methods</xref>
), we have arranged for the optimal estimate at time <italic>t</italic>
to be computable in closed form (see below). We emphasize that this computation does <italic>not</italic>
approximate the firing statistics of the Poisson observations as Gaussian; see the section <bold>The optimal filtering distribution</bold>
for details.</p>
<p>The task of <italic>learning</italic>
is to acquire the parameters that make it possible to carry out the estimation task. These parameters correspond to a dynamical system (e.g., the state transition matrix and the covariance of the state transition noise), and a model of how the states of this system give rise to observations. In our network, however, the parameters that are learned directly are the synaptic connections (weights) and the bias of each unit’s response function; the correspondence with the parameters of the dynamical system is not transparent (we explore it below).</p>
</sec>
<sec id="sec005"><title>The network</title>
<p>The central idea behind training our network, the “recurrent, exponential-family harmonium” (rEFH), is the choice of input data. In particular, each input vector consists of both the proprioceptive response to the current joint angle, and the activities of the hidden units at the previous time step. (For the case of a controlled dynamical system, the input vector also contains a noisy copy of the efferent motor command.) Biologically, this could be implemented via an additional population that simply reports, at a single time-step delay, the activities of the hidden units (<xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1B</xref>
, heavy black arrows; see also the section <bold>Cortical implementation</bold>
below). Conceptually, this choice of input data reflects the fact that filtering can be expressed as a “multisensory integration,” not between (e.g.) proprioceptive and visual inputs (cf. [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
]), but between proprioceptive inputs and a running best estimate of the state. We hypothesize that, because the hidden units learn to extract all the information available in their inputs [<xref rid="pcbi.1004554.ref009" ref-type="bibr">9</xref>
], the hidden vector at time <italic>t</italic>
− 1 will accumulate the information of the filtering distribution, <inline-formula id="pcbi.1004554.e003"><alternatives><graphic xlink:href="pcbi.1004554.e003.jpg" id="pcbi.1004554.e003g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M3"><mml:mrow><mml:mi>p</mml:mi>
<mml:mo>(</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold-italic">θ</mml:mi>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mn>0</mml:mn>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. Then at the next time step, the network will “integrate” this information with the current proprioceptive information about joint angle. (See <xref ref-type="supplementary-material" rid="pcbi.1004554.s002">S2 Text</xref>
for a longer discussion of this point.)</p>
<fig id="pcbi.1004554.g001" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g001</object-id>
<label>Fig 1</label>
<caption><title>Dynamical system and the neural network that learns it.</title>
<p>(A) The first four (of 1000) time steps of the linear system. The evolution is second-order, but only the current position (joint angle) is reported by the population of sensory neurons (orange)—fifteen of them, although only three are shown. (B) The input (“observed”) data for the harmonium (see text) at each time step are the current sensory activities, and recurrent activity of the hidden units, which amounts to a copy (heavy arrows) of their activity at the previous moment in time. (C) The joint angle. (D) (Top) Fifteen seconds of a typical trajectory (black) and the trajectory decoded from the position-encoding sensory population (orange), and from the hidden units (turquoise). (Bottom) The activity of the fifteen sensory neurons reported as a heat map. (E) Error statistics for angle decoding under different models (see text). Error bars mark the first and third quartile across twelve <italic>de novo</italic>
trained and tested networks of each type (see <xref ref-type="sec" rid="sec004">Methods</xref>
for training details).</p>
</caption>
<graphic xlink:href="pcbi.1004554.g001"></graphic>
</fig>
</sec>
<sec id="sec006"><title>Experiments</title>
<p>We therefore test our network by decoding, for all <italic>t</italic>
, hand position at time <italic>t</italic>
, from its hidden units. Rather than directly compare this estimate to the optimal estimate, which would provide no sense of scale, we compute error statistics for both. That is, we take the difference between the network’s estimate and the <italic>true</italic>
hand location; compute the mean and variance of this error across time steps (0 to <italic>T</italic>
= 1000) and trajectories (<italic>N</italic>
<sub>traj</sub>
= 40); and then compare these statistics for the rEFH and the optimum (OPT). We also compare error statistics from a “naïve” decoder (PROP) that simply decodes the current proprioceptive population, ignoring dynamics. It is the optimal decoder for data with no temporal dependencies.</p>
<p>The rEFH had to <italic>learn</italic>
to solve the estimation problem, and in practice, learned solutions will always be somewhat suboptimal, because of finite, noisy data and a nonconvex problem space. Therefore, a perhaps more useful point of comparison is the set of the error statistics from another model that has been trained on the same data. In particular, it is possible in certain cases to derive optimal parameter-update equations for a learning procedure (“expectation-maximization,” EM) that is guaranteed to reach at least local optima. Again by design of the dynamical system, and although the rEFH is not in theory limited to such data, such update equations are available (they are derived in <xref ref-type="supplementary-material" rid="pcbi.1004554.s001">S1 Text</xref>
). We therefore generate error statistics for this model (EM), as well. We emphasize, however, that EM was given additional information not provided to the rEFH: the order of the dynamical system, as well as the parameters of the observation model. The latter includes the best estimate of the stimulus given the population response, and the reliability of that “observation”—whereas the rEFH had to learn how to infer these values from the population response itself. Since EM is sensitive to the initial (random) values of the parameters it is to estimate, we present results for the best model from 20 random restarts (see <xref ref-type="sec" rid="sec004">Methods</xref>
); the same was done for the rEFH. To determine what order of dynamics the rEFH has learned, we also compare against lower-order models trained with EM. The order of these models is denoted with a superscript (e.g., EM<sup>2</sup>
).</p>
</sec>
<sec id="sec007"><title>Uncontrolled dynamical system</title>
<p>The generative model for the data appears in <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1A</xref>
. Conceptually, the problem is to track the shoulder joint (<xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1C</xref>
). To encourage second-order behavior, the parameters of the system were chosen make it underdamped (as e.g. when the arm hangs downward and acts as a pendulum; see <xref ref-type="sec" rid="sec004">Methods</xref>
), yielding trajectories like those shown in black in <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1D</xref>
, and the proprioceptive responses shown in orange. The rEFH was trained as a density estimator on these responses and the recurrent activity of its own hidden units at the previous time step (<xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1B</xref>
).</p>
<p>Error statistics for the various decoders are shown in <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1E</xref>
. Performance of the rEFH exceeds that of the naïve, purely sensory decoder (PROP), and approaches that of the optimum and the (second-order) EM-trained model (EM<sup>2</sup>
). The rEFH also outperforms the first-order, EM-trained model, EM<sup>1</sup>
, showing that it has learned to keep track not just of past positions, but of past velocities as well. We explore below <italic>how</italic>
it encodes position and velocity (<bold>Learned receptive fields and connectivity</bold>
).</p>
<p>We chose the dynamics of the underlying stimulus because they let us clearly see that the rEFH can learn a second-order model; that is, it learns to track the lawful changes in velocity, as well as position—even though only position information was available at each time step. But the results are robust across various dynamical models. <xref ref-type="fig" rid="pcbi.1004554.g002">Fig 2</xref>
shows results for 36 different dynamical systems which were created by varying the oscillator’s <xref ref-type="fig" rid="pcbi.1004554.g002">Fig 2A</xref>
stiffness, <xref ref-type="fig" rid="pcbi.1004554.g002">Fig 2B</xref>
damping, or <xref ref-type="fig" rid="pcbi.1004554.g002">Fig 2C</xref>
moment of inertia (colors as in <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1E</xref>
). For all of them, the rEFH outperforms the first-order model (EM<sup>1</sup>
), and performs close to the second-order model (EM<sup>2</sup>
).</p>
<fig id="pcbi.1004554.g002" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g002</object-id>
<label>Fig 2</label>
<caption><title>Mean squared errors (MSEs) for various dynamical systems.</title>
<p>To show the flexibility of the rEFH, we train one for each of several different second-order dynamical systems. Each subplot shows the results for twelve different systems in which a single parameter has been varied across a range of values; otherwise the systems are identical to the undriven model of <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1</xref>
. For each of those twelve dynamical systems, 20 rEFHs were trained (each with 150 hidden units), MSEs calculated, and the best selected (turquoise). The same was then done for EM<sup>1</sup>
(light red) and EM<sup>2</sup>
(blue), i.e., Kalman filters trained with EM, assuming either first- or second-order dynamics. The vertical black line in each plot indicates the dynamical system used in <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1</xref>
in the main text. (A) Varying the spring constant. The leftmost datum (<italic>k</italic>
= 0) corresponds to the “no-spring” model from which the RFs are analyzed below (although with lower transition noise); at this point, the dynamics can be well approximated by a first-order model. (B) Varying the damping coefficient. The spike in EM<sup>2</sup>
at <italic>c</italic>
= 0.4545 indicates that EM failed to find the second-order solution in any of its 20 attempts. As <italic>c</italic>
increases and the systems approach critical damping, however, first-order approximations are increasingly adequate. (C) Varying the moment of inertia.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g002"></graphic>
</fig>
</sec>
<sec id="sec008"><title>Controlled dynamical system</title>
<p>We now consider a system with inputs. Whereas the uncontrolled dynamical system, above, corresponds to the case where the arm is moved only passively by external forces, the controlled system corresponds to the more general case of self-motion. Controls are issued to the muscles of the arm by motor cortex, but a copy of these efferent signals is also fed back to posterior parietal cortex, <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3A</xref>
[<xref rid="pcbi.1004554.ref021" ref-type="bibr">21</xref>
]. This efference copy, being transmitted by a population of neurons, is assumed to be a noisy representation of the true control. Nor, presumably, is its role as a control signal explicitly given; rather, the network must learn, without supervision, to interpret it as such. This is precisely the learning problem faced by our network model (<xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1B</xref>
).</p>
<fig id="pcbi.1004554.g003" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g003</object-id>
<label>Fig 3</label>
<caption><title>Dynamical system with efference copy.</title>
<p>A reprise of <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1</xref>
for the controlled system. The evolution is now third-order, since the control signal is (close to) a random walk. (A) Angle (but not angular velocity) and the control signal are each noisily reported by populations of neurons (orange and purple; two of the fifteen neurons are depicted). (B) The training data for this harmonium at each time step are these two populations, and recurrent activity from the hidden units. (C) (Top) Fifteen seconds of a typical trajectory (black) and the trajectory decoded from the sensory population (orange), and from the hidden units (turquoise). The influence of the control can be seen in the (mildly) increasing amplitude of the oscillation. (Bottom) The activity of the fifteen sensory neurons is reported as a heat map. (D) Error statistics for joint-angle decoding under different models (see text). (E) The same as in (C), but the neurons carry “efference copy” rather than sensory information. (F) Error statistics for control-signal decoding under different models (see text). Error bars mark the first and third quartile across twelve <italic>de novo</italic>
trained and tested networks of each type.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g003"></graphic>
</fig>
<p>Realistic controls are correlated through time, so the control signal was given its own dynamics: a random walk with a very mild decay towards zero (see <xref ref-type="sec" rid="sec004">Methods</xref>
). This resulted in trajectories like that in <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3E</xref>
. The effect on angle can be seen in the corresponding trajectory of <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3C</xref>
: here, the control is driving the trajectory increasingly negative (cf. the first and second trough) in spite of the damping and the restoring force. In general, since the control is a random walk rather than merely white noise, the changes it effects on the position trajectory tend to accumulate.</p>
<p>Mean squared errors (MSEs) for the various filters of the controlled system are shown in <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3D</xref>
. The naïve model that ignores dynamics (PROP) is again the worst, as expected. Here the optimal EM-trained model is third-order (EM<sup>3</sup>
), since the second-order dynamical system is driven by a control with first-order dynamics. Again the rEFH performs close to this model, and outperforms the best lower-order model (EM<sup>2</sup>
). No trained model quite matches the true model’s performance (OPT), which result appears to be robust (see error bars), and presumably owes to the shape of the objective function (e.g., the optimal solution may be separated from suboptimal local maxima by deep valleys of low likelihood solutions).</p>
<p>So the rEFH has learned a third-order system. However, this does not <italic>per se</italic>
show that it has learned to use the efference-copy population; it might, for example, simply attribute all input to the system as white noise. To demonstrate that it does learn to use the controls, we compare it to the <italic>best</italic>
state estimates that can be made without efference copy. To produce such estimates, we fit a sixth filter, OBS, via regression, with full access to the state as well as the proprioceptive responses, but forced to assume zero control input. The resulting performance (<xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3D</xref>
, yellow bar) is clearly inferior to the other dynamical models.</p>
<p>Since the control signal is, like the state, only noisily observed by the network (via the efference copy), it is sensible to ask how well it can be decoded from the various models as well. And since the signal has its own (first-order) dynamics, it is possible for these models to make better estimates of the control at time <italic>t</italic>
than can be made from the efference copy alone at <italic>t</italic>
. <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3F</xref>
shows that this is indeed the case. All the filters perform about equally well, in comparison with the non-dynamical decoding of the efference-copy population (“EfCp”), although the harmonium is slightly inferior.</p>
</sec>
<sec id="sec009"><title>Distributions of performance across initializations and hidden-layer size</title>
<p>One known limitation on the learning capacity of the harmonium is the number of hidden units: the network requires sufficient representational power in this vector to encode the cumulants of the input distribution—in this case, the filter distribution. To determine what network size is necessary for maximal performance, we test a series of networks with systematically increasing numbers of hidden units. Because, however, the number of recurrent units must equal the number of hidden units, the ratio of hidden units to input units (recurrent, proprioceptive, and efference-copy) is necessarily upper bounded at unity. This asymptote presumably diminishes the returns of additional hidden units, even beyond those limitations imposed by the learning algorithm or the difficulty of the filtering task.</p>
<p>For each of twelve sizes, we train 20 networks <italic>de novo</italic>
, test them on a single fixed data set, and compute mean squared errors. The results are shown in <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4</xref>
, where network sizes (abscissae) are given by the number of hidden units. For the uncontrolled network, <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4A</xref>
, MSE of the best network (out of 20) clearly diminishes, asymptotically, with increasing numbers of hidden units. Beyond about 180 hidden units, no improvement in MSE is produced. For the controlled network, which has more to learn, the effect is even more severe (<xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4B</xref>
): increasing the number of hidden units beyond about 180 results in decreasing performance for even the best-performing networks at each network size. This suggests a limit to the complexity of the dynamical systems learned by this architecture, or with this learning procedure (for details of which, see <xref ref-type="sec" rid="sec004">Methods</xref>
).</p>
<fig id="pcbi.1004554.g004" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g004</object-id>
<label>Fig 4</label>
<caption><title>Box-and-whisker plot of MSEs for EM-based models and for rEFHs of various sizes.</title>
<p>Each box corresponds to 20 networks trained <italic>de novo</italic>
and tested on a common data set. Median MSE for each model is marked with a horizontal line; the box contains the interquartile range; whiskers extend to 1.5× the interquartile range, beyond which outliers are marked with plus signs. Performance among the EM-learned linear dynamical systems (the first two boxes, light red and blue) varies comparatively little, although large outliers are sometimes produced. In fact, the higher-order (best performing) models are all outliers. To facilitate comparison with the rEFHs, a line extends from the <italic>best</italic>
EM-based models across the entire plot. The remaining twelve boxes (turquoise) are rEFHs with different numbers of hidden units, listed on the abscissae. All have the same number of recurrent units as hidden units, and have a fixed number of “sensory units”: (A) 15 proprioceptive, or (B) 15 proprioceptive and 15 efference-copy. (A) The uncontrolled dynamical system. Overall, MSEs decline with increasing number of hidden/recurrent units, but appear to asymptote by 180 hidden units (ratio = 180/15 = 12). (B) The controlled dynamical system. The optimal number of hidden units is about 180 (ratio = 180/30 = 6), after which mean and variance (across networks) of MSE increases.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g004"></graphic>
</fig>
<p>On the other hand, rEFH learning appears to be more robust than EM learning, as can be seen in the box plots for the EM-based models (<xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4</xref>
). As with the rEFHs, 20 of each of the four EM-based models were trained from scratch, resulting in the distributions shown in light red and blue in <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4A and 4B</xref>
, (the narrowness of some of these distributions results in some very thin boxes). Although the EM algorithm guarantees convergence, it is only to a local (rather than global) optimum; which optimum is determined by the (random) initial parameters. The lower-order EM benchmarks (light red boxes; EM<sup>1</sup>
in <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4A</xref>
and EM<sup>2</sup>
in <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4B</xref>
) do indeed learn robustly, achieving nearly the same performance for all initializations. But the models with the true dynamical order (blue; EM<sup>2</sup>
in <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4A</xref>
and EM<sup>3</sup>
in <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4B</xref>
) exhibit a large performance distribution. These models are capable of outperforming the rEFH, but the runs that do are outliers. Thus, the large majority of the true-order EM-based models, for both the controlled and uncontrolled dynamical system, perform only about as well as their lower-order counterparts, which is inferior to all 20 of the 180-unit rEFH models. These rEFHs show comparatively little variation in performance—although that variance increases with the number of hidden units beyond 180.</p>
</sec>
</sec>
<sec id="sec010"><title>Learned receptive fields and connectivity</title>
<sec id="sec011"><title>How does the rEFH track the state?</title>
<p>Optimal (or nearly optimal) position estimation for these dynamical systems requires tracking velocity and position, so we plot receptive fields (RFs) in position-velocity space. Now, for oscillatory dynamics, high speeds rarely co-occur with positions far from zero (equilibrium), which leaves the “corners” of such RFs empty. This obscures the pattern of RFs and the corresponding state-estimation scheme learned by the rEFH. Therefore, for simplicity, we present results from a network trained on a third dynamical model (“no-spring”): uncontrolled, and with no spring force (see <xref ref-type="sec" rid="sec004">Methods</xref>
). (Similar results, albeit less clean, are observed in the corresponding analyses for oscillatory dynamics; see <xref ref-type="supplementary-material" rid="pcbi.1004554.s003">S3 Text</xref>
in the supporting material.) In <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5A</xref>
, the position-velocity receptive fields are plotted for all 225 hidden units of this rEFH, arranged in a 15 × 15 grid. The ordinate of each subsquare corresponds to position (increasing from top to bottom), and the abscissa to velocity (increasing from left to right). The large majority of receptive fields are negatively sloped “stripes” in this space. Interestingly, they resemble in this the receptive fields of neurons in MSTd of a rhesus macaque trained to track moving stimuli [<xref rid="pcbi.1004554.ref022" ref-type="bibr">22</xref>
]—although in that work there are positively-sloped stripes as well.</p>
<fig id="pcbi.1004554.g005" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g005</object-id>
<label>Fig 5</label>
<caption><title>Position and velocity receptive fields of hidden units.</title>
<p>In (A)-(C), pure white corresponds to a firing probability of one; pure black to zero. (A) Receptive fields for all 225 hidden units of the spring-free model (see text) in the space of (angular) position (ordinates) and velocity (abscissae). The angle limits and angular-velocity limits, indicated on the first (upper-left) receptive field, are the same for all units. (B) The predicted position-velocity receptive fields of units that have only the lagged-position tuning given by (C). The match with (A) is excellent for all but the anomalous 25 units at the right. (C) The same 225 units, each now plotted as a function of the time-lagged position with which that unit has maximal mutual information. Units have been arranged in order of increasing preferred position, whereas the units in (A) and (B) are arranged in order of maximally informative lags: from top to bottom and left to right, units are tuned for more temporally distant positions. This tuning gives rise to “stripes” in position-velocity space. For (A)-(C), the 25 units that do not appear to be well modeled by tuning to past positions have been placed at the end. (D) Histograms of the “preferred” lags, in terms both of time and (equivalently) discrete time steps, for five different networks. The normalized autocorrelation of the underlying dynamical system is superimposed. The central panel corresponds to the network analyzed in (A)-(C). The other four panels correspond to networks trained on observations from dynamical systems with different autocorrelations. From left to right panel, the dynamical systems get slower.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g005"></graphic>
</fig>
<p>Interpretation of these receptive fields is facilitated by an observation. If the velocity (<italic>ω</italic>
<sub><italic>t</italic>
</sub>
) is roughly constant over <italic>n</italic>
time steps of length Δ, then:
<disp-formula id="pcbi.1004554.e004"><alternatives><graphic xlink:href="pcbi.1004554.e004.jpg" id="pcbi.1004554.e004g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M4"><mml:mtable displaystyle="true"><mml:mtr><mml:mtd columnalign="right"><mml:msub><mml:mi>ω</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
</mml:mtd>
<mml:mtd columnalign="left"><mml:mrow><mml:mo>≈</mml:mo>
<mml:mfrac><mml:mrow><mml:msubsup><mml:mi>θ</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>−</mml:mo>
<mml:msubsup><mml:mi>θ</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mi>n</mml:mi>
<mml:mo>Δ</mml:mo>
</mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mrow><mml:mi>n</mml:mi>
<mml:mo>Δ</mml:mo>
</mml:mrow>
</mml:mfrac>
</mml:mrow>
</mml:mtd>
</mml:mtr>
<mml:mtr><mml:mtd columnalign="right"><mml:mrow><mml:mo>⇒</mml:mo>
<mml:msubsup><mml:mi>θ</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
</mml:mtd>
<mml:mtd columnalign="left"><mml:mrow><mml:mo>≈</mml:mo>
<mml:mi>n</mml:mi>
<mml:mo>Δ</mml:mo>
<mml:msub><mml:mi>ω</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>+</mml:mo>
<mml:msubsup><mml:mi>θ</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mi>n</mml:mi>
<mml:mo>Δ</mml:mo>
</mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:mtd>
</mml:mtr>
</mml:mtable>
</mml:math>
</alternatives>
</disp-formula>
the equation of a line in position-velocity space. Hence, such fields could be produced by neurons tuned simply for position at a delay, where the size of the delay (<italic>n</italic>
Δ) determines the slope (negative because we have plotted position as increasing from top to bottom, to match the corresponding figure in [<xref rid="pcbi.1004554.ref022" ref-type="bibr">22</xref>
]), and the “preferred” position determines the y-intercept. The equation is exact for <italic>n</italic>
= 1; but to the extent that velocity is not constant, the receptive fields will be diminished—as seen in the more irregular and faded character of receptive fields with greater slopes.</p>
<p>We therefore re-plot the receptive fields as a function of position only, but each at the time delay that maximizes mutual information between position and that hidden unit’s response (<xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5C</xref>
). Units are ordered by preferred position (whereas units in <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5A</xref>
were ordered by time delay). The resulting position tuning curves appear to tile space uniformly, with roughly constant receptive-field widths, suggesting that this is a concise description of the tuning curves that captures the computation being performed. As a final method of verification, we use these lagged-position receptive fields to generate idealized position-velocity receptive fields (<xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5B</xref>
; see <bold>Tuning analysis</bold>
for details.) The match with <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5A</xref>
is apparent.</p>
<p>This result is, perhaps, unexpected. It is possible to represent (an estimate of) the current state of a second-order dynamical system compactly by encoding current position and the current velocity. But it is also possible to represent it—seemingly less efficiently—in terms of the past positions alone, with the weight on each position decaying exponentially as a function of the number of time steps into the past. The rEFH appears to have learned a representation of this second type.</p>
<p>To determine if the weighting function applied to past positions by the rEFH does indeed correspond to such a scheme, we examine the distribution of “preferred” lags across hidden units (<xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5D</xref>
, center panel). Unsurprisingly, most units are tuned to the recent past, with an apparently monotonic decline into the past. Superimposed is the autocorrelation of the dynamical system on which the network was trained (see <xref ref-type="sec" rid="sec004">Methods</xref>
), normalized to have the same integral as the histogram. Evidently, the distribution of lags is well tuned to the dynamics. To confirm the robustness of this finding, we trained four new rEFHs on four different dynamical systems, identical to the one under discussion up to the damping coefficient (frictional force). From left to right in <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5D</xref>
, the dynamical systems were increasingly damped, resulting in longer autocorrelations (thick black lines). The temporal tuning of the rEFHs trained on these dynamical systems appears well matched in each case.</p>
</sec>
<sec id="sec012"><title>Responsiveness to instantaneous reliability</title>
<p>Thus, the rEFH tracks stimuli by encoding its position at various lags, with the number of units assigned to each lag decreasing exponentially with distance into the past, according to the autocorrelation of the signal. What is nevertheless not clear from <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5D</xref>
is whether the rEFH’s weighting of past positions takes into account the <italic>instantaneous</italic>
reliability of the proprioceptive encoding of joint angle. In our models, as (presumably) in the brain, instantaneous reliability varies because the proprioceptive report of joint angle is corrupted by Poisson noise. For Poisson spiking, the reliability—i.e., the inverse variance of the posterior distribution over joint location, conditioned on the spiking of all proprioceptive neurons—is proportional to the total number of spikes produced by the population at that time step. In the results discussed above, we also increased the fluctuations in this number by additionally varying the “gain” of the proprioceptive population, i.e., the single parameter that sets the height of all the tuning curves (see <bold>Input-data generation</bold>
). This is meant to model other random changes in the reliability of the proprioceptive report. Thus the optimal weighting of past position information, although on average an exponentially decaying function of time delay, will at any particular moment vary as a function of the recent, unpredictable, history of reliabilities: higher weights should be assigned to those time steps when the proprioceptive units had a collectively higher average firing rate, and <italic>vice versa</italic>
. A network that ignores such fluctuations will perform well, but suboptimally, perhaps explaining the (small) discrepancy between the MSEs for the rEFH and EM<sup>2</sup>
seen in <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1E</xref>
.</p>
<p>Here we show that the rEFH is indeed sensitive to instantaneous reliability. We retest the network of the section <bold>Uncontrolled dynamical system</bold>
on <italic>noiseless</italic>
sensory data, i.e., using the mean spike counts of the proprioceptive neurons rather than Poisson samples drawn from those means. This allows us to set the total spike count essentially directly. In this noiseless test, a higher total spike count does not correspond to a more reliable signal: the signal is perfectly reliable for all spike counts. Hence for any total spike count, minimal (zero) error could be achieved by a “filter” that relied entirely on the current sensory input. Nevertheless, the optimal filter <italic>for the data on which the rEFH was trained</italic>
, viz., OPT, will <italic>not</italic>
rely entirely upon this current sensory information, but rather will weight it in proportion to the total spike count. In consequence, OPT will achieve lower MSEs on data with greater total spike counts, since for such data it will lean more heavily on the perfect sensory information. This is the pattern observed in <xref ref-type="fig" rid="pcbi.1004554.g006">Fig 6A</xref>
. If the rEFH is also treating total spike count correctly—i.e., if it is properly sensitive to instantaneous reliability—its MSEs will exhibit a similar pattern. And indeed, they do (<xref ref-type="fig" rid="pcbi.1004554.g006">Fig 6B</xref>
).</p>
<fig id="pcbi.1004554.g006" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g006</object-id>
<label>Fig 6</label>
<caption><title>Network sensitivity to instantaneous reliability.</title>
<p>The instantaneous (one-time-step) reliability of sensory information is determined by the total number of spikes across the sensory population within one time step. An optimal filter will up-weight sensory information that is more reliable (and vice versa). If such a filter is run on <italic>noiseless</italic>
sensory data, then its errors will be smaller for sensory input with more total spikes (higher gain), since it will up-weight the perfect sensory information. (A) Box-and-whisker plot (interpretation as in <xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4</xref>
) of mean squared errors for the optimal model (OPT), when tested on noiseless sensory data and a range of gains. For each gain on the abscissa, the filter was tested on 12 sets of 320 trajectories apiece, for which the sensory gain was fixed throughout. Higher-gain trajectories yield lower mean errors, as expected. (B) The same plot for the network (rEFH). The magnitude of the MSEs is larger than for the optimal filter, as in <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1E</xref>
, but the pattern is the same, showing that the rEFH has indeed learned to treat higher-gain (more-spikes) sensory information as more reliable.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g006"></graphic>
</fig>
<p>Here we emphasize again that, for the optimal filter (as well as OBS and EM<sup><italic>n</italic>
</sup>
), we provided the transformation from total spike count to sensory reliability, whereas the the rEFH had to learn this transformation. Likewise, in engineered solutions to tracking problems, the Kalman filter is usually simplified to learn a single, average reliability for all time. We have demonstrated that the rEFH is not similarly limited, since it can learn to use instantaneous indicators of reliability if they are present in the observations.</p>
</sec>
<sec id="sec013"><title>Emergence of receptive fields</title>
<p>Since velocity is not reported directly by the sensory (“proprioceptive”) population, the network will not immediately develop tuning for it. <xref ref-type="fig" rid="pcbi.1004554.g007">Fig 7</xref>
illustrates its emergence across training trials. (We return here to the “no-spring” model for comparison with <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5A</xref>
.) Since position information is a useful “feature” for explaining the proprioceptive inputs, the hidden units learn to extract it after just 100 batches of training (<xref ref-type="fig" rid="pcbi.1004554.g007">Fig 7A</xref>
). But information that is in the hidden units will also appear in the input, at a one step time delay, via the recurrent units (see <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1B</xref>
). So at this point in training, the input contains information about both the current position (in the proprioceptive population) and about the past position (in the recurrent population). This makes extraction of velocity information possible. It is <italic>useful</italic>
because the stimuli obey second-order dynamics: knowing the relationship between past and present position allows each to provide information about the other, yielding overall superior estimates. And indeed, by 200 batches of learning, some velocity tuning appears, evidenced by the sloping of receptive fields in position-velocity space, <xref ref-type="fig" rid="pcbi.1004554.g007">Fig 7B</xref>
. But these units look back no more than a few steps in time. By 1000 batches (<xref ref-type="fig" rid="pcbi.1004554.g007">Fig 7C</xref>
), strong velocity tuning is evident, although the full distribution of lags (slopes) has yet to emerge (cf. <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5A</xref>
, which is after some 100,000 batches of training).</p>
<fig id="pcbi.1004554.g007" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g007</object-id>
<label>Fig 7</label>
<caption><title>Emergence of position and velocity receptive during training.</title>
<p>Velocity tuning takes longer to emerge than position tuning, because velocity information is only available in the inputs after position has already been learned by the hidden units—and subsequently fed back. (A) After 100 batches of training, stripes are mostly horizontal or very shallowly sloped in position-velocity space—no velocity tuning. (B) By 200 batches, velocity tuning is evident across most units. (C) Later, at 1000 batches, the slopes of the “stripes” have increased, indicating position tuning for more temporally distant (past) stimuli.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g007"></graphic>
</fig>
</sec>
<sec id="sec014"><title>Organization of the learned weight matrix</title>
<p>The network learns to model the dynamics by making changes to the synaptic connection strengths, summarized by the weight matrices <italic>W</italic>
<sub>prop</sub>
(sensory to hidden) and <italic>W</italic>
<sub>fb</sub>
(recurrent to hidden). To understand better the mechanism of the network we examine these matrices (<xref ref-type="fig" rid="pcbi.1004554.g008">Fig 8</xref>
), again for the “no-spring” model. (The corresponding figure for the non-zero stiffness model, slightly more difficult to interpret, is shown in <xref ref-type="supplementary-material" rid="pcbi.1004554.s003">S3 Text</xref>
.) In the arbitrarily ordered form in which they are learned, they are difficult to decipher, but interesting patterns emerge when they are reordered by the parameters of the receptive fields. The reordering is applied to both the hidden and the recurrent units. (The original, topographical ordering of the sensory units is retained.)</p>
<fig id="pcbi.1004554.g008" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g008</object-id>
<label>Fig 8</label>
<caption><title>The weight matrices.</title>
<p>The organization of connections between inputs (sensory or recurrent units) and the hidden layer can be visualized by sorting units by (A) preferred stimulus angle or (B) by preferred lag (increasing from lower left to upper right). Here we analyze the “no-spring” network (see text). In both subfigures, both the recurrent and hidden units have been re-sorted; the sensory units remain organized by preferred angle. Note that self-connections (along the diagonal) are in fact more than 0.5, but the plot saturates at this value to make the other connections more visible.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g008"></graphic>
</fig>
<p>First, reordering by the hidden units’ preferred stimulus angles (“PA”), <xref ref-type="fig" rid="pcbi.1004554.g008">Fig 8A</xref>
, reveals that the difference in PA between two hidden units dictates the sign of their connection. Hidden units with similar PAs have positive connections, and units with “out of phase” (recall that the stimulus is a circular variable) PAs have negative connections. This results from the continuity of the stimulus trajectories: the network “expects” the stimulus to move from any given position (encoded in the recurrent units) to a nearby one (encoded in the hidden units).</p>
<p>Second, we reorder the units by “preferred lag,” <italic>τ</italic>
—i.e., the time delay at which mutual information between sensory input and hidden-unit activity is maximal (<xref ref-type="fig" rid="pcbi.1004554.g008">Fig 8B</xref>
). Again, units with similar <italic>τ</italic>
are preferentially wired together.</p>
</sec>
</sec>
<sec id="sec015"><title>Cortical implementation</title>
<p>More than one cortical area is thought to subserve object tracking. Since we have in this study focused on the task of tracking one’s own limbs, we consider posterior parietal cortex (PPC), which is thought to be responsible for this task [<xref rid="pcbi.1004554.ref010" ref-type="bibr">10</xref>
, <xref rid="pcbi.1004554.ref011" ref-type="bibr">11</xref>
]. The computation may well be distributed across the PPC, but we focus on just one that has been particularly implicated [<xref rid="pcbi.1004554.ref011" ref-type="bibr">11</xref>
], Brodmann Area 5. Our aim is to show that our neural network and its learning scheme are consistent with what is known about the connectivity of Area 5, both interlaminar and inter-areal. In particular, we consider its connections with the primary motor area (M1) and primary somatosensory cortex (S1). Our proposed implementation is speculative and not the only one possible; e.g., we identify the “recurrent” units with another layer of Area 5, but they might alternatively correspond to another area of PPC.</p>
<p><xref ref-type="fig" rid="pcbi.1004554.g009">Fig 9A</xref>
summarizes the training procedure from an algorithmic perspective (see <xref ref-type="sec" rid="sec004">Methods</xref>
for details). In <xref ref-type="fig" rid="pcbi.1004554.g009">Fig 9B</xref>
, as in <xref ref-type="fig" rid="pcbi.1004554.g009">Fig 9A</xref>
, input comes from two sources. Feedforward, proprioceptive input (<inline-formula id="pcbi.1004554.e005"><alternatives><graphic xlink:href="pcbi.1004554.e005.jpg" id="pcbi.1004554.e005g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M5"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) from primary somatosensory cortex, S1 (especially BA3a), projects to layer IV [<xref rid="pcbi.1004554.ref023" ref-type="bibr">23</xref>
]. A copy of the efferent command (<inline-formula id="pcbi.1004554.e006"><alternatives><graphic xlink:href="pcbi.1004554.e006.jpg" id="pcbi.1004554.e006g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M6"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) feeds back from M1 to layer I of Area 5 [<xref rid="pcbi.1004554.ref023" ref-type="bibr">23</xref>
]. Layer II/III of Area 5 in turn projects forward to M1 [<xref rid="pcbi.1004554.ref024" ref-type="bibr">24</xref>
]. Layer I is not believed to contain cell bodies [<xref rid="pcbi.1004554.ref025" ref-type="bibr">25</xref>
], so we take these to be the terminal branches of the apical dendrites of layer II/III cells (which are also lightly labeled by anterograde tracers injected in M1 [<xref rid="pcbi.1004554.ref023" ref-type="bibr">23</xref>
]). Within Area 5, we propose that the temporally delayed recurrency (<bold>Z</bold>
<sub><italic>t</italic>
−1</sub>
) of the rEFH is provided by the loop from layer II/III down to VI, then up to V, before modulating the activity of layer II/III neurons, consistent with the anatomy of Area 5 [<xref rid="pcbi.1004554.ref025" ref-type="bibr">25</xref>
]. Layer IV and III, as well as V and III, also have reciprocal connections [<xref rid="pcbi.1004554.ref025" ref-type="bibr">25</xref>
], as required for the rEFH training procedure. The latter loop has in fact been hypothesized to give rise to rhythmic activity in rat parietal cortex [<xref rid="pcbi.1004554.ref026" ref-type="bibr">26</xref>
].</p>
<fig id="pcbi.1004554.g009" orientation="portrait" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1004554.g009</object-id>
<label>Fig 9</label>
<caption><title>Training and testing: in the model, and in a cortical implementation.</title>
<p>(A) The training and testing procedure in the model. Three discrete times steps are arrayed vertically. At each one, the current arm position (<bold>θ</bold>
<sub><italic>t</italic>
</sub>
) is reported by the proprioceptive population (orange) with Poisson-distributed spike counts, as shown in the first column. The current motor command is likewise reported by an “efference-copy” population (purple). Second column: this spiking, along with recurrent activities (dark turqoise), stochastically drives single spikes in the hidden layer (turquoise). Third column: these spikes in turn drive the three “input” populations (this is not required during testing); a “copy” of the hidden vector is also saved to serve as recurrent activity at the next time step (curved black arrow). Fourth column: finally, the input populations drive the hidden layer once more, after which the weights are changed according to <xref ref-type="disp-formula" rid="pcbi.1004554.e040">Eq. 9</xref>
(also not required during testing). At every time step, the current joint angle and current control are decoded naïvely from the current activities of their respective input populations; with Kalman filters that are recursively updated based on these activities (not depicted in this figure); and from the hidden units of the rEFH. (B) A possible cortical implementation of the rEFH. The cortical layers of Brodmann Area 5 and its inputs are identified with elements of the rEFH by color. The synaptic connections are denoted by the arrows and their corresponding weight matrices. Primary somatosensory cortex (S1) provides feedforward proprioceptive input to layer IV, while primary motor cortex (M1) provides feedback input—a copy of the efferent command—to layer II/III. The recurrent signal of the rEFH (heavy curved arrow in (A)) is identified with a reverberatory loop from II/III to VI, to V, and back to II/III.</p>
</caption>
<graphic xlink:href="pcbi.1004554.g009"></graphic>
</fig>
<p>According to the learning and filtering schemes of our model, the temporal flow of information is as follows. Sensory input (<inline-formula id="pcbi.1004554.e007"><alternatives><graphic xlink:href="pcbi.1004554.e007.jpg" id="pcbi.1004554.e007g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M7"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) and efference copy (<inline-formula id="pcbi.1004554.e008"><alternatives><graphic xlink:href="pcbi.1004554.e008.jpg" id="pcbi.1004554.e008g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M8"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) arrive at, respectively, layer IV of BA5 and the feedback layer (presumably VI) of M1. At the same time, a “copy” (which could be any information-preserving transformation) of activity from layer II/III of BA5 (<bold>z</bold>
<sub><italic>t</italic>
-1</sub>
) passes down to layer V. Next, the spiking in these layers (M1 layer VI, BA5 layer IV, BA5 layer V) drives spiking (<bold>z</bold>
<sub><italic>t</italic>
</sub>
) in BA5 layer II/III. These responses encode, according to the model, the optimal estimate of the limb, and this information will ultimately become the temporally delayed recurrent activities identified above. For learning, however, it is also necessary that this activity drive spiking in M1 (<inline-formula id="pcbi.1004554.e009"><alternatives><graphic xlink:href="pcbi.1004554.e009.jpg" id="pcbi.1004554.e009g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M9"><mml:msubsup><mml:mrow><mml:mover accent="true"><mml:mi mathvariant="bold">r</mml:mi>
<mml:mo>^</mml:mo>
</mml:mover>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
), BA5 layer IV (<inline-formula id="pcbi.1004554.e010"><alternatives><graphic xlink:href="pcbi.1004554.e010.jpg" id="pcbi.1004554.e010g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M10"><mml:msubsup><mml:mrow><mml:mover accent="true"><mml:mi mathvariant="bold">r</mml:mi>
<mml:mo>^</mml:mo>
</mml:mover>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
), and BA5 layer V (<inline-formula id="pcbi.1004554.e011"><alternatives><graphic xlink:href="pcbi.1004554.e011.jpg" id="pcbi.1004554.e011g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M11"><mml:msub><mml:mover accent="true"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mi>t</mml:mi>
</mml:msub>
</mml:math>
</alternatives>
</inline-formula>
), through the reciprocal connectivity lately noted. A “copy” of the layer II/III activity (<bold>z</bold>
<sub><italic>t</italic>
</sub>
) is simultaneously propagated down to layer VI. Lastly, the activities in M1, BA5 layer IV, and BA5 layer V again drive activity (<inline-formula id="pcbi.1004554.e012"><alternatives><graphic xlink:href="pcbi.1004554.e012.jpg" id="pcbi.1004554.e012g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M12"><mml:msub><mml:mover accent="true"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mi>t</mml:mi>
</mml:msub>
</mml:math>
</alternatives>
</inline-formula>
) in BA5 layer II/III. At the same time, the “copy” of layer II/III activity (<bold>z</bold>
<sub><italic>t</italic>
</sub>
) is communicated up to layer V.</p>
</sec>
</sec>
<sec sec-type="conclusions" id="sec016"><title>Discussion</title>
<sec id="sec017"><title>Summary of results</title>
<p>We have shown that a neural network (the “rEFH”) with a biologically plausible architecture and synaptic-plasticity rule can learn to track moving stimuli, in the sense that its downstream (“hidden”) units learn to encode (nearly) the most accurate estimate of stimulus location possible, given the entire history of incoming sensory information (Figs <xref ref-type="fig" rid="pcbi.1004554.g001">1</xref>
and <xref ref-type="fig" rid="pcbi.1004554.g002">2</xref>
). This requires learning a model of the stimulus dynamics. This is (as far as we know) the first biologically plausible model that has been shown to learn to solve this task. Moreover, the network learns the <italic>reliability</italic>
of the sensory signal: the trained network leans more heavily on the internal model when the sensory signal is less reliable, and more heavily on the sensory signal when it is more reliable (<xref ref-type="fig" rid="pcbi.1004554.g006">Fig 6</xref>
).</p>
<p>We are particularly interested in tracking the state of one’s own limbs. Here, additional information about stimulus location is thought to be available in the form of a “copy,” relayed to the posterior parietal cortex, of the efferent motor command [<xref rid="pcbi.1004554.ref021" ref-type="bibr">21</xref>
]. And indeed, when such signals are available to our network, it learns to make use of them appropriately to track the arm more precisely—in spite of the fact that none of the incoming signals is “labeled” according to its role (<xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3</xref>
). Although an expectation-maximization (EM) algorithm can sometimes learn a Kalman filter that noticeably outperforms the best rEFH on these data, it usually does not (<xref ref-type="fig" rid="pcbi.1004554.g004">Fig 4</xref>
). That is, learning in the rEFH is more robust than EM in the sense that the variance in performance across models trained <italic>de novo</italic>
is smaller, albeit at the price of a bias towards worse models. Finally, and surprisingly, the downstream neurons of the trained network track a moving stimulus by encoding its position at various time lags (<xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5</xref>
).</p>
</sec>
<sec id="sec018"><title>Related work</title>
<p>The earliest implementation of dynamical state estimation (“filtering”) in neural architecture comes from Rao and Ballard [<xref rid="pcbi.1004554.ref027" ref-type="bibr">27</xref>
]. Their model, like ours, assigns a central role to recurrent connections, but as predictive coders rather than simply delayed copies of previous neural states. Likewise, the network connectivity is acquired with an unsupervised and local learning rule, a variant on EM. However, the authors do not train their network on moving objects or moving images, presumably because convergence of the neural state under their learning scheme is slow compared with any plausible stimulus dynamics. Instead, the connectivity is acquired on static images. Performance on state-estimation tasks is not tested.</p>
<p>Several groups have hand wired neural networks to act as state estimators [<xref rid="pcbi.1004554.ref007" ref-type="bibr">7</xref>
, <xref rid="pcbi.1004554.ref008" ref-type="bibr">8</xref>
, <xref rid="pcbi.1004554.ref028" ref-type="bibr">28</xref>
]. Although these papers do not address our central concern, the learning problem, it is nevertheless useful to compare the resulting architectures with our rEFH. For example, Beck and colleagues constructed a neural network to implement the Kalman filter on linear probabilistic population codes, as in this work, and showed its performance (measured in information loss) to be nearly optimal. From analytical considerations, the authors showed that the required operations on neural firing rates are weighted summation (as in our network) and a quadratic operation (that acts like a divisive normalization in the steady state). In our rEFH, on the other hand, the only nonlinearities are elementwise: interaction between inputs is always in the form of a weighted sum. That the rEFH can nevertheless filter (nearly) optimally is possible because we do not require, as they do, that the output population encode information in the same way as the inputs (sc., that the posterior distribution over the stimulus have linear sufficient statistics; see <xref ref-type="supplementary-material" rid="pcbi.1004554.s006">S4 Text</xref>
). This critical difference provides the basis for an experimental discrimination between the respective models. Likewise, filters have been hand wired into attractor networks [<xref rid="pcbi.1004554.ref028" ref-type="bibr">28</xref>
] and spike-based (rather than rate-based) networks [<xref rid="pcbi.1004554.ref008" ref-type="bibr">8</xref>
]. The latter in particular argues that the precise arrival time of spikes contains information about the stimulus, rather than the average rate across time, as in in our model.</p>
<p>An approach that does include learning comes from Huys, Zemel, Natarajan, and Dayan [<xref rid="pcbi.1004554.ref029" ref-type="bibr">29</xref>
, <xref rid="pcbi.1004554.ref030" ref-type="bibr">30</xref>
]. The authors formulate the problem in terms very similar to ours, but they allow more general dynamical systems generated by Gaussian processes, and the basic unit of information is spikes rather than spike counts (although approximations that ignore precise arrival times lose little information [<xref rid="pcbi.1004554.ref029" ref-type="bibr">29</xref>
]). The most significant difference with our work is that the authors learn the parameters of their network with a supervised, non-local rule, which they do not consider to be a biological mechanism. But again the comparison is instructive. We are able to formulate an unsupervised rule because we approach the filtering problem indirectly: Natarajan and colleagues require the posterior distribution, conditioned on hidden-unit activities, to be factorizable over hidden-unit spikes (so that a third layer can consider those spikes separately), and then force it to match the true filtering distribution by directly descending the KL divergence between them [<xref rid="pcbi.1004554.ref030" ref-type="bibr">30</xref>
]. We, on the other hand, force the network to be a good model of its incoming data—which, when some of those data are past hidden-unit activities, achieves the same end.</p>
<p>In the machine-learning literature, Hinton and colleagues have proposed three variants on a theme quite similar to ours [<xref rid="pcbi.1004554.ref031" ref-type="bibr">31</xref>
–<xref rid="pcbi.1004554.ref033" ref-type="bibr">33</xref>
], although different in important ways. Most importantly, in all three, the past hidden-unit activities are treated by the learning rule as (fixed) biases rather than as input data; i.e., they cannot be modified during the “down pass” of contrastive-divergence training. That these activities ought to be treated as data, we argue more rigorously in a forthcoming work.</p>
<p>The earliest variant [<xref rid="pcbi.1004554.ref031" ref-type="bibr">31</xref>
], the “spiking Boltzmann machine,” is, like ours, a temporal extension of the restricted Boltzmann machine that is trained with the contrastive-divergence rule. Hidden units are directly influenced by past hidden-unit activities, as with the rEFH, but possibly from temporal distances <italic>τ</italic>
that are greater than one time step (<italic>contra</italic>
the rEFH). However, the weights from a particular “past” hidden unit at various delays (e.g., from <inline-formula id="pcbi.1004554.e013"><alternatives><graphic xlink:href="pcbi.1004554.e013.jpg" id="pcbi.1004554.e013g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M13"><mml:mrow><mml:msubsup><mml:mi>z</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mi>n</mml:mi>
<mml:mi>τ</mml:mi>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mi>n</mml:mi>
<mml:mo>∈</mml:mo>
<mml:mrow><mml:mo>{</mml:mo>
<mml:mn>1</mml:mn>
<mml:mo>,</mml:mo>
<mml:mn>2</mml:mn>
<mml:mo>,</mml:mo>
<mml:mo>.</mml:mo>
<mml:mo>.</mml:mo>
<mml:mo>.</mml:mo>
<mml:mo>}</mml:mo>
</mml:mrow>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
to <inline-formula id="pcbi.1004554.e014"><alternatives><graphic xlink:href="pcbi.1004554.e014.jpg" id="pcbi.1004554.e014g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M14"><mml:msubsup><mml:mi>z</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>j</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) are constrained to be identical up to a fixed (not learned) exponential decay. The motivation was to model the influence of past spikes in a biologically plausible way: Whereas in our rEFH, the (one-time-step delayed) past hidden activities are maintained in a separate population of neurons (<xref ref-type="fig" rid="pcbi.1004554.g009">Fig 9</xref>
), in the “spiking Boltzmann machine” their effect on current hidden units is interpreted simply as the decaying influence of their original arrival. This makes it plausible, unlike in the rEFH, to include influences at delays greater than one time step. On the other hand, it necessitates treating those effects as biases rather than data.</p>
<p>It is difficult to judge the limitations this imposes on the model, since the authors do not quantify its performance. However, they do investigate more thoroughly performance of a similar, but more powerful network. The “temporal restricted Boltzmann machine” (TRBM) [<xref rid="pcbi.1004554.ref032" ref-type="bibr">32</xref>
] is a spiking Boltzmann machine without the constraint that the weights decay exponentially backwards in time; instead, they are learned freely and independently for all time. The order of the dynamical system that can be learned by this network turns out, unlike ours, to be tied to <italic>τ</italic>
: TRBMs with <italic>τ</italic>
= 1 (like the rEFH) can learn only random-walk behavior (first-order dynamics) [<xref rid="pcbi.1004554.ref033" ref-type="bibr">33</xref>
]. This can (presumably) be overcome by including connections back as many time steps as the order of the system to be learned, but it is not obvious what biological mechanism could maintain copies of past activities at distant lags, or determine <italic>a priori</italic>
how many such lags to maintain. The same authors show that this problem can be alleviated with a variant architecture, the “recurrent temporal RBM” (RTRBM) [<xref rid="pcbi.1004554.ref033" ref-type="bibr">33</xref>
], but it requires a non-causal learning rule (backpropagation through time), again making it a poor model for neural function. For neither model do the authors precisely quantify its filtering performance; we do in a forthcoming study.</p>
</sec>
<sec id="sec019"><title>Implications of the model</title>
<p>Our simulations demonstrate three things: First, the rEFH is capable of learning to “track” moving stimuli, i.e. to estimate their dynamical state, and nearly as well as an optimal algorithm, as has been seen behaviorally in humans [<xref rid="pcbi.1004554.ref034" ref-type="bibr">34</xref>
]. In fact, the network learns to encode the full posterior distribution over the stimulus, rather than just its peak: although we did not show it directly, it must, since the variance of this (Gaussian) distribution is required to combine properly the previous best estimate with the current sensory information. And rather than relying on a fixed estimate of sensory reliability, the network learns to take into account instantaneous changes in it (<xref ref-type="fig" rid="pcbi.1004554.g006">Fig 6</xref>
).</p>
<p>Second, the network does not require a special architecture or <italic>ad hoc</italic>
modifications. It is, rather, identical, up to the choice of input populations, to the network and learning rule in our previous work [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
]. Thus, if the input populations are proprioceptive and <italic>recurrent</italic>
units, it will learn to estimate dynamical state; if they also include efference copy, it will learn the influence of motor commands on stimulus dynamics. If they are proprioceptive and <italic>visual</italic>
reports of a common stimulus, it will learn to perform multisensory integration; if a gaze-angle-reporting population is also present, to transform the visual signal by that angle before integrating (“coordinate transformations”); if the stimulus distribution is non-uniform, to encode that distribution [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
]. (We have shown elsewhere, in terms of information theory, why this is the case [<xref rid="pcbi.1004554.ref009" ref-type="bibr">9</xref>
]. For further discussion of the relationship between the static and dynamical computations, see <xref ref-type="supplementary-material" rid="pcbi.1004554.s002">S2 Text</xref>
.) Thus, the network provides a very general model for posterior parietal cortex, where some combination of all of these signals is often present.</p>
<p>Third, the model makes some predictions about the encoding scheme, receptive fields, and connectivity of cortical areas that track objects. As with all models, we take certain elements of ours to be essential and others to be adventitious. That learning in posterior-parietal circuits can be well described as a form of latent-variable density estimation, for example, is central to our theory; but the precise form of the learning rule (“one-step contrastive divergence”), although plausible, is not. Our theory requires that sensory neurons encode distributions over stimulus position, but the representation scheme need not be probabilistic population codes of the Pougetian variety [<xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
]. Here we list three predictions that do follow from essential aspects of the network.
<list list-type="order"><list-item><p>The network learns to track by encoding past positions. This is a non-obvious scheme (it is not, e.g., the one used by the Kalman filter) and apparently results from the fact that only position information is reported by the sensory afferents. It is possible that such receptive fields (<xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5A</xref>
) are in fact found in MSTd of monkeys that have been trained to track moving objects [<xref rid="pcbi.1004554.ref022" ref-type="bibr">22</xref>
]. Now, in many circuits, velocity is detected at early stages. But even when velocity is directly reported by the inputs to an rEFH, tuning to past positions still appears, albeit with lower prevalence (see <xref ref-type="supplementary-material" rid="pcbi.1004554.s003">S3 Text</xref>
).</p>
<p>More generally, we predict that higher derivatives (e.g., acceleration), especially those not directly available in sensory input, will be encoded via delayed, lower derivatives (e.g., velocities)—as long as those higher-order states have lawful dynamics.</p>
</list-item>
<list-item><p>During learning, receptive fields for position emerge before those for velocity. This is a necessary consequence of density estimation on recurrent units. A similar proviso attaches: where velocity information is directly reported, it is acceleration-coding that will emerge over time.</p>
</list-item>
<list-item><p>The use of delayed, feedback connections in neural circuits is a mechanism for learning dynamical properties of stimuli. Under this prediction, primary sensory areas that process information with very little temporal structure—e.g., smell—will lack the dense feedback found in, e.g., visual areas. Alternatively, the recurrency might be identified with interlaminar, rather than interareal, structure, as we have hypothesized (<xref ref-type="fig" rid="pcbi.1004554.g009">Fig 9B</xref>
)—which would explain why piriform cortex only needs three layers.</p>
</list-item>
</list>
</p>
</sec>
<sec id="sec020"><title>Neural computation in posterior parietal cortex</title>
<p>More generally, our investigation was motivated by two main ideas. The first is that populations of neurons, in virtue of their natural variability, encode probability distributions over stimuli (rather than point estimates) [<xref rid="pcbi.1004554.ref001" ref-type="bibr">1</xref>
, <xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
]. Encoding certainty or “reliability” is a necessity for optimal integration of dynamic sensory information, since it determines the relative weight given to (a) current sensory information and (b) the prediction of the internal model. But rather than <italic>explicitly</italic>
encoding the reliability of the stimulus location—e.g., via neurons that are “tuned” to reliability, as other neurons are tuned to location itself—this reliability is identified with the inverse variance of the posterior distribution over the stimulus, <italic>p</italic>
(<bold>s</bold>
∣<bold>r</bold>
), conditioned on the population activity [<xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
]. This distribution arises as a natural consequence of the (putative) fact that neural responses are noisy, and can therefore be characterized by a likelihood, <italic>p</italic>
(<bold>r</bold>
∣<bold>s</bold>
) [<xref rid="pcbi.1004554.ref001" ref-type="bibr">1</xref>
]. If reliability were not encoded this way, our learning scheme would not work: it would have no way of knowing what to do with those reliabilities, which would be to it indistinguishable from (e.g.) the location of another stimulus.</p>
<p>The second idea is that higher sensory areas, like posterior parietal cortex and MSTd, can encode more precise distributions over the location (e.g.) of a stimulus than that provided by their sensory afferents at any given moment in time. This is due, essentially, to the continuity of the physical world: at successive moments in time, objects tend to remain near their previous positions. More precise localizations can consequently be achieved by a form of averaging that, because objects do move, accounts for the predictable changes in position from moment to moment. This requires learning a model of those predictable changes. The rich statistical structure of the sensory afferents—including efference copy of motor commands that may be influencing the evolution of the stimulus to be tracked, as when tracking one’s own limbs—makes it possible to learn the model from those inputs alone. This unsupervised learning is a much more efficient approach than trying to use the few bits of information that may be available in the form of reward: very few rewards can be reaped before an animal can control its own limbs.</p>
<p>In the special case of linear dynamics and Gaussian noise, these two problems—learning a dynamical model, and filtering in that model—have known algorithmic solutions: an expectation-maximization algorithm and the Kalman filter, respectively. Rather than try to map operations on vectors and matrices directly onto neural activity and learning rules, we have taken a more general approach, showing how a rather general neural-network architecture that tries to build good models for its inputs can learn to solve the problem, if those inputs are suitably chosen: temporally delayed recurrent activity from downstream units must be among the inputs. Our network learns by a local, Hebbian rule operating on spike-count correlations, although it remains to relate these to more specific biological learning rules, like STDP.</p>
</sec>
</sec>
<sec sec-type="materials|methods" id="sec021"><title>Methods</title>
<p>Notation is standard: capital letters for random variables, lowercase for their realizations; boldfaced font for vectors, italic for scalars. Capitalized italics are also used for matrices (context distinguishes them from random scalars).</p>
<sec id="sec022"><title>Input-data generation</title>
<p>We describe the most general dynamical system and observation model to be learned: a controlled, second-order, discrete-time, stochastic, linear dynamical system, whose “observations” or outputs come in the form of linear probabilistic population codes [<xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
]; cf. <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3A</xref>
. The uncontrolled model of the section <bold>Uncontrolled dynamical system</bold>
(<xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1A</xref>
) is a special case (see below). We interpret the plant to be a rotational joint, so distance is in units of radians; and the control to be a torque, hence in Joules/radian.</p>
<p>The primary rationale for our choice of dynamics and observation model was to show what kinds of computational issues the recurrent, exponential-family harmonium (rEFH) can overcome—issues which it must overcome if it is to be a good model for the way cortex learns to solves the problem. In particular, it might appear that the rEFH can learn relationships only between its current inputs and the previous ones, since its recurrent inputs are from the previous time step only (see <xref ref-type="fig" rid="pcbi.1004554.g009">Fig 9A</xref>
). Therefore, we let the inputs report position only, but make the (hidden) dynamics second-order: velocity, as well as position, depends on previous position and velocity. If the rEFH can learn to associate only current and previous inputs, it can learn only first-order dynamics from these data. Furthermore, to clearly distinguish models that have learned second-order dynamics from those that have learned only a first-order approximation, we let the true dynamics be a (damped) oscillator (first-order systems cannot oscillate). Although the demonstration is in terms of positions and velocities, the point is more general: if the rEFH can learn second-order dynamics from position reports, it can learn higher temporal dynamics from lower-order data more generally.</p>
<p>The controlled, single-joint limb obeys:
<disp-formula id="pcbi.1004554.e015"><alternatives><graphic xlink:href="pcbi.1004554.e015.jpg" id="pcbi.1004554.e015g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M15"><mml:mrow><mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">θ</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>+</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">θ</mml:mi>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mi mathvariant="script">N</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:mi>A</mml:mi>
<mml:msub><mml:mi mathvariant="bold-italic">θ</mml:mi>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>+</mml:mo>
<mml:mi mathvariant="bold">b</mml:mi>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">μ</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msub><mml:mo>Σ</mml:mo>
<mml:mi>θ</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(1)</label>
</disp-formula>
where the vector random variable <bold>Θ</bold>
<sub><italic>t</italic>
</sub>
consists of angle and angular velocity. The control signal (torque) has itself first-order dynamics:
<disp-formula id="pcbi.1004554.e016"><alternatives><graphic xlink:href="pcbi.1004554.e016.jpg" id="pcbi.1004554.e016g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M16"><mml:mrow><mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>+</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mi mathvariant="script">N</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:mi>α</mml:mi>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi>μ</mml:mi>
<mml:mi>u</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mi>σ</mml:mi>
<mml:mi>u</mml:mi>
<mml:mn>2</mml:mn>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(2)</label>
</disp-formula>
making the combined system third-order. The initial state and control are also normally distributed:
<disp-formula id="pcbi.1004554.e017"><alternatives><graphic xlink:href="pcbi.1004554.e017.jpg" id="pcbi.1004554.e017g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M17"><mml:mrow><mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">θ</mml:mi>
<mml:mn>0</mml:mn>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mn>0</mml:mn>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mi mathvariant="script">N</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msub><mml:mi>ν</mml:mi>
<mml:mn>0</mml:mn>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msub><mml:mi>ϒ</mml:mi>
<mml:mn>0</mml:mn>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>.</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(3)</label>
</disp-formula>
The current (time <italic>t</italic>
) joint position and control are noisily encoded in the spike counts of populations of neurons, whose Gaussian-shaped tuning curves (<italic>f</italic>
<sub><italic>i</italic>
</sub>
) smoothly tile their respective spaces, proprioceptive (angle) and control (torque). Spike counts are drawn from (conditionally) independent Poisson distributions:
<disp-formula id="pcbi.1004554.e018"><alternatives><graphic xlink:href="pcbi.1004554.e018.jpg" id="pcbi.1004554.e018g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M18"><mml:mrow><mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msub><mml:mi mathvariant="bolditalic">θ</mml:mi>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mrow><mml:msubsup><mml:mi>g</mml:mi>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:munder><mml:mo>∏</mml:mo>
<mml:mi>i</mml:mi>
</mml:munder>
<mml:mtext>Pois</mml:mtext>
<mml:mrow><mml:mo>[</mml:mo>
<mml:mrow><mml:msubsup><mml:mi>r</mml:mi>
<mml:mrow><mml:mi>i</mml:mi>
<mml:mo>,</mml:mo>
<mml:mi>t</mml:mi>
</mml:mrow>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mrow><mml:mo>|</mml:mo>
</mml:mrow>
<mml:msubsup><mml:mrow><mml:msubsup><mml:mi>g</mml:mi>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:msub><mml:mi>f</mml:mi>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mrow><mml:mo>(</mml:mo>
<mml:mi>C</mml:mi>
<mml:msub><mml:mi mathvariant="bolditalic">θ</mml:mi>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:mrow>
<mml:mo>]</mml:mo>
</mml:mrow>
<mml:mo>,</mml:mo>
<mml:mspace width="18.06749pt"></mml:mspace>
<mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mrow><mml:msubsup><mml:mi>g</mml:mi>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:munder><mml:mo>∏</mml:mo>
<mml:mi>i</mml:mi>
</mml:munder>
<mml:mtext>Pois</mml:mtext>
<mml:mrow><mml:mo>[</mml:mo>
<mml:mrow><mml:msubsup><mml:mi>r</mml:mi>
<mml:mrow><mml:mi>i</mml:mi>
<mml:mo>,</mml:mo>
<mml:mi>t</mml:mi>
</mml:mrow>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mrow><mml:mo>|</mml:mo>
</mml:mrow>
<mml:msubsup><mml:mrow><mml:msubsup><mml:mi>g</mml:mi>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mspace width="1pt"></mml:mspace>
<mml:msub><mml:mi>f</mml:mi>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mrow><mml:mo>(</mml:mo>
<mml:mi>h</mml:mi>
<mml:msubsup><mml:mi>u</mml:mi>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:mrow>
<mml:mo>]</mml:mo>
</mml:mrow>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(4)</label>
</disp-formula>
with <italic>C</italic>
= [1 0] and <italic>h</italic>
= 1. Here the <italic>g</italic>
<sub><italic>t</italic>
</sub>
are “gains,” scaling factors for the mean spike count [<xref rid="pcbi.1004554.ref002" ref-type="bibr">2</xref>
, <xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
]. Because the signal-to-noise ratio increases with mean for Poisson random variables, these gains essentially scale (linearly) the reliability of each population. Therefore, in order to model instant-to-instant changes in sensory reliability, the gains of each population were chosen independently and uniformly:
<disp-formula id="pcbi.1004554.e019"><alternatives><graphic xlink:href="pcbi.1004554.e019.jpg" id="pcbi.1004554.e019g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M19"><mml:mrow><mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mrow><mml:msubsup><mml:mi>g</mml:mi>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mi mathvariant="script">U</mml:mi>
<mml:mo>(</mml:mo>
<mml:mn>6</mml:mn>
<mml:mo>.</mml:mo>
<mml:mn>4</mml:mn>
<mml:mo>,</mml:mo>
<mml:mn>9</mml:mn>
<mml:mo>.</mml:mo>
<mml:mn>6</mml:mn>
<mml:mo>)</mml:mo>
<mml:mo>,</mml:mo>
<mml:mspace width="36.135pt"></mml:mspace>
<mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mrow><mml:msubsup><mml:mi>g</mml:mi>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mi mathvariant="script">U</mml:mi>
<mml:mo>(</mml:mo>
<mml:mn>6</mml:mn>
<mml:mo>.</mml:mo>
<mml:mn>4</mml:mn>
<mml:mo>,</mml:mo>
<mml:mn>9</mml:mn>
<mml:mo>.</mml:mo>
<mml:mn>6</mml:mn>
<mml:mo>)</mml:mo>
<mml:mo>.</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(5)</label>
</disp-formula>
Since the discrete time interval for a single draw from <xref ref-type="disp-formula" rid="pcbi.1004554.e018">Eq. 4</xref>
is 0.05 s (see below), these gains correspond to <italic>maximal</italic>
firing rates between 130 and 192 spikes/second, reasonable rates for neurons in cortex. The joint distribution of the states, controls, their observations, and the gains is the product of Eqs <xref ref-type="disp-formula" rid="pcbi.1004554.e015">1</xref>
–<xref ref-type="disp-formula" rid="pcbi.1004554.e019">5</xref>
, multiplied across all time.</p>
<p>In accordance with the broad tuning of higher sensory areas, the “standard deviation,” <italic>σ</italic>
<sub>tc</sub>
, of the tuning curves,
<disp-formula id="pcbi.1004554.e020"><alternatives><graphic xlink:href="pcbi.1004554.e020.jpg" id="pcbi.1004554.e020g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M20"><mml:mrow><mml:msub><mml:mi>f</mml:mi>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mrow><mml:mo>(</mml:mo>
<mml:mi>x</mml:mi>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mtext>exp</mml:mtext>
<mml:mo>{</mml:mo>
<mml:mo>−</mml:mo>
<mml:mfrac><mml:msup><mml:mrow><mml:mo>(</mml:mo>
<mml:mi>x</mml:mi>
<mml:mo>−</mml:mo>
<mml:msub><mml:mi>ξ</mml:mi>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mn>2</mml:mn>
</mml:msup>
<mml:mrow><mml:mn>2</mml:mn>
<mml:msubsup><mml:mi>σ</mml:mi>
<mml:mtext>tc</mml:mtext>
<mml:mn>2</mml:mn>
</mml:msubsup>
</mml:mrow>
</mml:mfrac>
<mml:mo>}</mml:mo>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</disp-formula>
was chosen so that the full-width at half maximum is one-sixth of the space of feasible joint angles/torques, for all preferred stimuli <italic>ξ</italic>
<sub><italic>i</italic>
</sub>
. However, joints and torques can in fact leave these “feasible spaces”: Although the system was designed to be stable (eigenvalues of the state-transition matrix are within the unit circle), trajectories are nevertheless unbounded, since the input noise is unbounded (normally distributed). We chose not to impose hard joint and torque limits, because this would make the dynamics nonlinear, vitiating the optimality calculations. Instead, stimuli beyond the feasible space simply “wrap” onto the opposite side of encoding space; that is, each population tiles its corresponding stimulus <italic>modulo</italic>
the length of its feasible space.</p>
<p>But for the dynamical systems on which model performance was tested, parameters were chosen to make wrapping unlikely (but cf. the “no-spring” model described below). In particular, we used the discrete-time approximation to a damped harmonic oscillator, i.e., <inline-formula id="pcbi.1004554.e021"><alternatives><graphic xlink:href="pcbi.1004554.e021.jpg" id="pcbi.1004554.e021g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M21"><mml:mrow><mml:mi>m</mml:mi>
<mml:mover accent="true"><mml:mi>θ</mml:mi>
<mml:mo>¨</mml:mo>
</mml:mover>
<mml:mo>+</mml:mo>
<mml:mi>c</mml:mi>
<mml:mover accent="true"><mml:mi>θ</mml:mi>
<mml:mo>˙</mml:mo>
</mml:mover>
<mml:mo>+</mml:mo>
<mml:mi>k</mml:mi>
<mml:mi>θ</mml:mi>
<mml:mo>=</mml:mo>
<mml:mi>u</mml:mi>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
:
<disp-formula id="pcbi.1004554.e022"><alternatives><graphic xlink:href="pcbi.1004554.e022.jpg" id="pcbi.1004554.e022g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M22"><mml:mrow><mml:mi>A</mml:mi>
<mml:mo>=</mml:mo>
<mml:mo>[</mml:mo>
<mml:mtable><mml:mtr><mml:mtd><mml:mn>1</mml:mn>
</mml:mtd>
<mml:mtd><mml:mo>Δ</mml:mo>
</mml:mtd>
</mml:mtr>
<mml:mtr><mml:mtd><mml:mrow><mml:mo>−</mml:mo>
<mml:mfrac><mml:mi>k</mml:mi>
<mml:mi>m</mml:mi>
</mml:mfrac>
<mml:mo>Δ</mml:mo>
</mml:mrow>
</mml:mtd>
<mml:mtd><mml:mrow><mml:mn>1</mml:mn>
<mml:mo>−</mml:mo>
<mml:mfrac><mml:mi>c</mml:mi>
<mml:mi>m</mml:mi>
</mml:mfrac>
<mml:mo>Δ</mml:mo>
</mml:mrow>
</mml:mtd>
</mml:mtr>
</mml:mtable>
<mml:mo>]</mml:mo>
<mml:mo>,</mml:mo>
<mml:mspace width="18.06749pt"></mml:mspace>
<mml:mi mathvariant="bold">b</mml:mi>
<mml:mo>=</mml:mo>
<mml:mo>[</mml:mo>
<mml:mtable><mml:mtr><mml:mtd><mml:mn>0</mml:mn>
</mml:mtd>
</mml:mtr>
<mml:mtr><mml:mtd><mml:mfrac><mml:mo>Δ</mml:mo>
<mml:mi>m</mml:mi>
</mml:mfrac>
</mml:mtd>
</mml:mtr>
</mml:mtable>
<mml:mo>]</mml:mo>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</disp-formula>
with moment of inertia <italic>m</italic>
= 5 J⋅ s<sup>2</sup>
/rad<sup>2</sup>
, viscous damping <italic>c</italic>
= 0.25 J⋅ s/rad<sup>2</sup>
, ideal-spring stiffness <italic>k</italic>
= 3 J/rad<sup>2</sup>
, and sampling interval Δ = 0.05 s. This makes the system stable and underdamped (oscillatory). The control decay, <italic>α</italic>
, in <xref ref-type="disp-formula" rid="pcbi.1004554.e016">Eq 2</xref>
was set to 0.9994, making the dynamics close to a random walk, but mildly decaying towards zero.</p>
<p>These parameters and the noise variances were chosen so that the system could not be well approximated by a lower-order one—i.e., so that the uncontrolled and controlled systems were “truly” second- and third-order (respectively). This was accomplished by ensuring that the Hankel singular values [<xref rid="pcbi.1004554.ref014" ref-type="bibr">14</xref>
] for the system, with output matrix <italic>C</italic>
= [1 0] and input matrix set by the noise variances, were within one order of magnitude of each other; that is, ensuring that the transfer function from noise to joint angle had roughly equal power in all modes. For the uncontrolled system, this was achieved with Σ<sub><italic>θ</italic>
</sub>
= diag([5<sc>e</sc>
-7, 5<sc>e</sc>
-5]); for the controlled system, Σ<sub><italic>θ</italic>
</sub>
= diag([5<sc>e</sc>
-5, 1<sc>e</sc>
-6]) and <inline-formula id="pcbi.1004554.e023"><alternatives><graphic xlink:href="pcbi.1004554.e023.jpg" id="pcbi.1004554.e023g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M23"><mml:mrow><mml:msubsup><mml:mi>σ</mml:mi>
<mml:mi>u</mml:mi>
<mml:mn>2</mml:mn>
</mml:msubsup>
<mml:mo>=</mml:mo>
<mml:mn>7</mml:mn>
<mml:mo>.</mml:mo>
<mml:mn>5</mml:mn>
<mml:mtext>E</mml:mtext>
<mml:mo>−</mml:mo>
<mml:mn>4</mml:mn>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. While this last choice of noise is large enough to ensure that the control’s contribution to the dynamics is significant, it is also small enough to keep wrapping rare. This facilitates the comparison between the benchmark models (see below), which are acquired from non-wrapped trajectories, and the rEFH, which learns from sensory inputs with periodic tuning curves. That is, for fast enough trajectories on a circle, the dynamics would no longer be locally linear, and the learning and filtering tasks no longer comparable.</p>
<p>The only other difference between the uncontrolled and controlled dynamical systems was that the former had, of course, no control signal (or simply <bold>b</bold>
= <bold>0</bold>
) and no control observations (efference copy). For all models, the bias terms were set to zero: <bold><italic>μ</italic>
</bold>
<sub><italic>θ</italic>
</sub>
= <bold>0</bold>
and <italic>μ</italic>
<sub><italic>u</italic>
</sub>
= 0.</p>
<p>The initial positions for all trajectories were drawn from a uniform distribution across joint space (shoulder <italic>θ</italic>
∈ [−<italic>π</italic>
/3, <italic>π</italic>
/3] radians; <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1C</xref>
), up to a margin of 0.05 radians from the joint limits (to discourage state transitions out of the feasible space); for EM learning (see below), this was treated as an infinite-covariance Gaussian centered in the middle of joint space. The initial velocity and initial control were normally distributed very tightly about zero, with a standard deviation of <inline-formula id="pcbi.1004554.e024"><alternatives><graphic xlink:href="pcbi.1004554.e024.jpg" id="pcbi.1004554.e024g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M24"><mml:mrow><mml:msqrt><mml:mn>5</mml:mn>
</mml:msqrt>
<mml:mtext>E</mml:mtext>
<mml:mo>−</mml:mo>
<mml:mn>5</mml:mn>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
(rad/s and J/rad, resp.). Hence <bold><italic>ν</italic>
</bold>
<sub>0</sub>
= <bold>0</bold>
, <bold>Υ</bold>
<sub>0</sub>
= diag([∞, 5<sc>e</sc>
-10, 5<sc>e</sc>
-10]). The range (modulus) of feasible controls is <italic>u</italic>
∈ [−1.25, 1.25] J/rad.</p>
<p>For the receptive-field (RF) analyses, we used a third dynamical system. In the harmonic oscillator, whether driven or undriven, the non-zero stiffness (<italic>k</italic>
above) couples velocity to position, making high speeds and far-from-zero positions unlikely to co-occur. This makes the RF analysis unreliable in the “corners” of position-velocity space, and the overall velocity-encoding harder to interpret. For the analyses presented in Figs <xref ref-type="fig" rid="pcbi.1004554.g005">5</xref>
, <xref ref-type="fig" rid="pcbi.1004554.g007">7</xref>
and <xref ref-type="fig" rid="pcbi.1004554.g008">8</xref>
, therefore, we trained a (third) rEFH on a simplified version (“no-spring”) of the uncontrolled dynamics, setting the spring constant to zero (eliminating oscillations). To encourage full exploration of the space, the variance of the state-transition noise was also increased by a factor of 50. The more and less autocorrelated variants of <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5D</xref>
were created by simply scaling up or down the damping coefficient: from left to right, <italic>c</italic>
= 0.25/4, 0.25/2, 0.25, 0.25 * 2, 0.25 * 4. For completeness, we nevertheless include, in the Supplement, the harder-to-interpret RF analyses for the rEFH trained on the (undriven) harmonic oscillator (<xref ref-type="supplementary-material" rid="pcbi.1004554.s003">S3 Text</xref>
).</p>
</sec>
<sec id="sec023"><title>The recurrent, exponential-family harmonium (rEFH)</title>
<p>The network is very similar to that in [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
], but we repeat the description here briefly. The harmonium is a generalization of the restricted Boltzmann machine (RBM) beyond Bernoulli units to other random variables in the exponential family [<xref rid="pcbi.1004554.ref003" ref-type="bibr">3</xref>
]. That is, it is a two-layer network with full interlayer connections and no intralayer connections, which can be thought of as a Markov random field (undirected graphical model) or as a neural network. In our implementation (see Figs <xref ref-type="fig" rid="pcbi.1004554.g001">1B</xref>
and <xref ref-type="fig" rid="pcbi.1004554.g003">3B</xref>
), hidden units (turquoise, <bold>Z</bold>
<sub><italic>t</italic>
</sub>
) and recurrent units (dark turqoise, <bold>Z</bold>
<sub><italic>t</italic>
−1</sub>
) are binary (spike/no spike), and the “proprioceptive” (orange, <inline-formula id="pcbi.1004554.e025"><alternatives><graphic xlink:href="pcbi.1004554.e025.jpg" id="pcbi.1004554.e025g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M25"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) and “efference-copy” (purple, <inline-formula id="pcbi.1004554.e026"><alternatives><graphic xlink:href="pcbi.1004554.e026.jpg" id="pcbi.1004554.e026g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M26"><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
) populations are non-negative integers (spike counts). For all networks, the number of recurrent units is the same as the number of downstream or “hidden” units, because recurrent units at time <italic>t</italic>
carry the activities of the hidden units at time <italic>t</italic>
− 1—making the harmonium recurrent through time (rEFH). We chose <italic>N</italic>
<sub>hid</sub>
= <italic>N</italic>
<sub>recurrent</sub>
= 240 for the network trained on the uncontrolled system, and <italic>N</italic>
<sub>hid</sub>
= <italic>N</italic>
<sub>recurrent</sub>
= 180 for the controlled system. We used fifteen proprioceptive units (<italic>N</italic>
<sub>prop</sub>
) and, for the network trained on the controlled system, fifteen efference-copy units (<italic>N</italic>
<sub>efcp</sub>
), so the total number of “observed” (or “input”) variables was 255 = <italic>N</italic>
<sub>recurrent</sub>
+ <italic>N</italic>
<sub>prop</sub>
for the uncontrolled model and 210 = <italic>N</italic>
<sub>recurrent</sub>
+ <italic>N</italic>
<sub>prop</sub>
+ <italic>N</italic>
<sub>efcp</sub>
for the controlled model.</p>
<p>During training and testing, the layers of the rEFH reciprocally drive each other, yielding samples from the following distributions:
<disp-formula id="pcbi.1004554.e027"><alternatives><graphic xlink:href="pcbi.1004554.e027.jpg" id="pcbi.1004554.e027g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M27"><mml:mrow><mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">Z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>∼</mml:mo>
<mml:mi>q</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:munderover><mml:mo>∏</mml:mo>
<mml:mi>i</mml:mi>
<mml:msub><mml:mi>N</mml:mi>
<mml:mtext>hid</mml:mtext>
</mml:msub>
</mml:munderover>
<mml:mtext>Bern</mml:mtext>
<mml:mo>[</mml:mo>
<mml:mrow><mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:mi>σ</mml:mi>
<mml:mo>(</mml:mo>
<mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msub><mml:mi>W</mml:mi>
<mml:mtext>fb</mml:mtext>
</mml:msub>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi>W</mml:mi>
<mml:mtext>prop</mml:mtext>
</mml:msub>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">W</mml:mi>
<mml:mtext>ctrl</mml:mtext>
</mml:msub>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi mathvariant="bold">b</mml:mi>
<mml:mtext>hid</mml:mtext>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>]</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(6a)</label>
</disp-formula>
<disp-formula id="pcbi.1004554.e028"><alternatives><graphic xlink:href="pcbi.1004554.e028.jpg" id="pcbi.1004554.e028g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M28"><mml:mrow><mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">Z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>∼</mml:mo>
<mml:mi>q</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:munderover><mml:mo>∏</mml:mo>
<mml:mi>i</mml:mi>
<mml:msub><mml:mi>N</mml:mi>
<mml:mtext>hid</mml:mtext>
</mml:msub>
</mml:munderover>
<mml:mtext>Bern</mml:mtext>
<mml:mo>[</mml:mo>
<mml:mrow><mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:mi>σ</mml:mi>
<mml:mo>(</mml:mo>
<mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msubsup><mml:mi>W</mml:mi>
<mml:mtext>fb</mml:mtext>
<mml:mtext>T</mml:mtext>
</mml:msubsup>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi mathvariant="bold">b</mml:mi>
<mml:mtext>fb</mml:mtext>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>]</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(6b)</label>
</disp-formula>
<disp-formula id="pcbi.1004554.e029"><alternatives><graphic xlink:href="pcbi.1004554.e029.jpg" id="pcbi.1004554.e029g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M29"><mml:mrow><mml:msubsup><mml:mi mathvariant="bold">R</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>∼</mml:mo>
<mml:mi>q</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:munderover><mml:mo>∏</mml:mo>
<mml:mi>i</mml:mi>
<mml:msub><mml:mi>N</mml:mi>
<mml:mtext>prop</mml:mtext>
</mml:msub>
</mml:munderover>
<mml:mtext>Pois</mml:mtext>
<mml:mo>[</mml:mo>
<mml:mrow><mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:mtext>exp</mml:mtext>
<mml:mo>(</mml:mo>
<mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msubsup><mml:mi>W</mml:mi>
<mml:mtext>prop</mml:mtext>
<mml:mtext>T</mml:mtext>
</mml:msubsup>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi mathvariant="bold">b</mml:mi>
<mml:mtext>prop</mml:mtext>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>]</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(6c)</label>
</disp-formula>
<disp-formula id="pcbi.1004554.e030"><alternatives><graphic xlink:href="pcbi.1004554.e030.jpg" id="pcbi.1004554.e030g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M30"><mml:mrow><mml:msubsup><mml:mi mathvariant="bold">R</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>∼</mml:mo>
<mml:mi>q</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:munderover><mml:mo>∏</mml:mo>
<mml:mi>i</mml:mi>
<mml:msub><mml:mi>N</mml:mi>
<mml:mtext>efcp</mml:mtext>
</mml:msub>
</mml:munderover>
<mml:mtext>Pois</mml:mtext>
<mml:mo>[</mml:mo>
<mml:mrow><mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msubsup><mml:mi mathvariant="bold">r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:mtext>exp</mml:mtext>
<mml:mo>(</mml:mo>
<mml:msub><mml:mrow><mml:mo>{</mml:mo>
<mml:msubsup><mml:mi>W</mml:mi>
<mml:mtext>efcp</mml:mtext>
<mml:mtext>T</mml:mtext>
</mml:msubsup>
<mml:msub><mml:mrow><mml:msubsup><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi mathvariant="bold">b</mml:mi>
<mml:mtext>efcp</mml:mtext>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>]</mml:mo>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(6d)</label>
</disp-formula>
which corresponds to Gibbs sampling from the joint distribution represented by the harmonium, <inline-formula id="pcbi.1004554.e031"><alternatives><graphic xlink:href="pcbi.1004554.e031.jpg" id="pcbi.1004554.e031g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M31"><mml:mrow><mml:mi>q</mml:mi>
<mml:mo>(</mml:mo>
<mml:msub><mml:mi mathvariant="bold">z</mml:mi>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msub><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>;</mml:mo>
<mml:mi>W</mml:mi>
<mml:mo>,</mml:mo>
<mml:mi mathvariant="bold">b</mml:mi>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. The letter <italic>q</italic>
is used for the probability density function assigned by the rEFH to distinguish it from the true distribution over the observed variables, <inline-formula id="pcbi.1004554.e032"><alternatives><graphic xlink:href="pcbi.1004554.e032.jpg" id="pcbi.1004554.e032g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M32"><mml:mrow><mml:mi>p</mml:mi>
<mml:mo>(</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">r</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. Here the notation {<bold>x</bold>
}<sub><italic>i</italic>
</sub>
means the <italic>i</italic>
<sup>th</sup>
element of the vector <bold>x</bold>
; the matrices <italic>W</italic>
and vectors <bold>b</bold>
are the synaptic connection strengths (“weights”) and biases, respectively; and the neural nonlinearities, the logistic (<italic>σ</italic>
(<italic>x</italic>
) = 1/(1 + <italic>e</italic>
<sup>−<italic>x</italic>
</sup>
)) and exponential funtions, were chosen to produce means for each distribution that are in the appropriate interval ([0, 1] and <inline-formula id="pcbi.1004554.e033"><alternatives><graphic xlink:href="pcbi.1004554.e033.jpg" id="pcbi.1004554.e033g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M33"><mml:mrow><mml:msub><mml:mi>ℝ</mml:mi>
<mml:mo>+</mml:mo>
</mml:msub>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
, resp.). The entire procedure is depicted graphically in <xref ref-type="fig" rid="pcbi.1004554.g009">Fig 9A</xref>
.</p>
<sec id="sec024"><title>Training</title>
<p>Although the ultimate goal of training is to make the network able to solve the filtering problem, this is achieved indirectly by making the harmonium a good model for the data on which it was trained. That is, the harmonium should assign probability (<italic>q</italic>
) to the observed data (<bold>Y</bold>
) that matches the probability with which they actually appear (<italic>p</italic>
); in short, the goal is to achieve:
<disp-formula id="pcbi.1004554.e034"><alternatives><graphic xlink:href="pcbi.1004554.e034.jpg" id="pcbi.1004554.e034g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M34"><mml:mrow><mml:mi>q</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">y</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>;</mml:mo>
<mml:mi>W</mml:mi>
<mml:mo>,</mml:mo>
<mml:mi mathvariant="bold">b</mml:mi>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">y</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(7)</label>
</disp-formula>
equality between the “model distribution” and “data distribution,” by adjustment of the weights and biases of the network. In our case, <bold>Y</bold>
= [<bold>Y</bold>
<sub>0</sub>
, ⋯, <bold>Y</bold>
<sub><italic>T</italic>
</sub>
], a collection of observations across time, where intuitively the observations at time <italic>t</italic>
are the responses of the proprioceptive and efference-copy populations, <inline-formula id="pcbi.1004554.e035"><alternatives><graphic xlink:href="pcbi.1004554.e035.jpg" id="pcbi.1004554.e035g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M35"><mml:mrow><mml:msub><mml:mrow><mml:mi mathvariant="bold">Y</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>=</mml:mo>
<mml:mrow><mml:mo>[</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">R</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>]</mml:mo>
</mml:mrow>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. However, these random variables are not independent across time; that is, <inline-formula id="pcbi.1004554.e036"><alternatives><graphic xlink:href="pcbi.1004554.e036.jpg" id="pcbi.1004554.e036g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M36"><mml:mrow><mml:mi>p</mml:mi>
<mml:mo stretchy="false">(</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>T</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>T</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo stretchy="false">)</mml:mo>
<mml:mo>≠</mml:mo>
<mml:mstyle displaystyle="true"><mml:msub><mml:mo>∏</mml:mo>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mi>p</mml:mi>
</mml:mstyle>
<mml:mo stretchy="false">(</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo stretchy="false">)</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. In order, then, to make possible incremental training—weight changes without first collecting population responses for all time, [0, …, <italic>T</italic>
]—we train on the augmented observation vector:
<disp-formula id="pcbi.1004554.e037"><alternatives><graphic xlink:href="pcbi.1004554.e037.jpg" id="pcbi.1004554.e037g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M37"><mml:mrow><mml:msub><mml:mi mathvariant="bold">Y</mml:mi>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>=</mml:mo>
<mml:mo stretchy="false">[</mml:mo>
<mml:msub><mml:mi mathvariant="bold">Z</mml:mi>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo stretchy="false">]</mml:mo>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(8)</label>
</disp-formula>
where <bold>Z</bold>
<sub><italic>t</italic>
−1</sub>
are the hidden-unit activities at the previous time step. Intuitively, the addition of these recurrent activities renders the data independent because they recursively accumulate all the information contained in their inputs [<xref rid="pcbi.1004554.ref009" ref-type="bibr">9</xref>
].</p>
<p>Weight changes are made proportional to the approximate gradient of a function (“one-step contrastive divergence,” CD<sub>1</sub>
) that has <xref ref-type="disp-formula" rid="pcbi.1004554.e034">Eq 7</xref>
at its minimum [<xref rid="pcbi.1004554.ref012" ref-type="bibr">12</xref>
, <xref rid="pcbi.1004554.ref013" ref-type="bibr">13</xref>
]. In exponential family harmoniums, following this gradient is particularly simple: Stimuli in the world drive the input populations (<bold>y</bold>
, <xref ref-type="disp-formula" rid="pcbi.1004554.e018">Eq 4</xref>
), which drive the hidden units (<bold>z</bold>
, <xref ref-type="disp-formula" rid="pcbi.1004554.e027">Eq 6a</xref>
), which reciprocally drive the input populations (<inline-formula id="pcbi.1004554.e038"><alternatives><graphic xlink:href="pcbi.1004554.e038.jpg" id="pcbi.1004554.e038g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M38"><mml:mover accent="true"><mml:mi>y</mml:mi>
<mml:mo>^</mml:mo>
</mml:mover>
</mml:math>
</alternatives>
</inline-formula>
, Eqs <xref ref-type="disp-formula" rid="pcbi.1004554.e028">6b</xref>
, <xref ref-type="disp-formula" rid="pcbi.1004554.e029">6c</xref>
and <xref ref-type="disp-formula" rid="pcbi.1004554.e030">6d</xref>
), which drive the hidden units once more (<inline-formula id="pcbi.1004554.e039"><alternatives><graphic xlink:href="pcbi.1004554.e039.jpg" id="pcbi.1004554.e039g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M39"><mml:mover accent="true"><mml:mi>z</mml:mi>
<mml:mo>^</mml:mo>
</mml:mover>
</mml:math>
</alternatives>
</inline-formula>
, <xref ref-type="disp-formula" rid="pcbi.1004554.e027">Eq 6a</xref>
); after which parameters are changed according to:
<disp-formula id="pcbi.1004554.e040"><alternatives><graphic xlink:href="pcbi.1004554.e040.jpg" id="pcbi.1004554.e040g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M40"><mml:mrow><mml:mo>Δ</mml:mo>
<mml:mi>W</mml:mi>
<mml:mo>∝</mml:mo>
<mml:mrow><mml:mi mathvariant="bold">y</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
<mml:msup><mml:mrow><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
</mml:mrow>
<mml:mtext>T</mml:mtext>
</mml:msup>
<mml:mo>−</mml:mo>
<mml:mover accent="true"><mml:mrow><mml:mi mathvariant="bold">y</mml:mi>
</mml:mrow>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:msup><mml:mover accent="true"><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mtext>T</mml:mtext>
</mml:msup>
<mml:mo>,</mml:mo>
<mml:mspace width="36.135pt"></mml:mspace>
<mml:mo>Δ</mml:mo>
<mml:msub><mml:mi mathvariant="bold">b</mml:mi>
<mml:mi>y</mml:mi>
</mml:msub>
<mml:mo>∝</mml:mo>
<mml:mrow><mml:mi mathvariant="bold">y</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
<mml:mo>−</mml:mo>
<mml:mover accent="true"><mml:mrow><mml:mi mathvariant="bold">y</mml:mi>
</mml:mrow>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mo>,</mml:mo>
<mml:mspace width="36.135pt"></mml:mspace>
<mml:mo>Δ</mml:mo>
<mml:msub><mml:mi mathvariant="bold">b</mml:mi>
<mml:mi>z</mml:mi>
</mml:msub>
<mml:mo>∝</mml:mo>
<mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mrow></mml:mrow>
<mml:mrow></mml:mrow>
<mml:mo>−</mml:mo>
<mml:mover accent="true"><mml:mrow><mml:mi mathvariant="bold">z</mml:mi>
</mml:mrow>
<mml:mo>^</mml:mo>
</mml:mover>
<mml:mo>.</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(9)</label>
</disp-formula>
Note that the learning rule is local and Hebbian (correlational). The entire procedure amounts to taking a full step of Gibbs sampling in a Markov chain that has been initialized at a vector sampled from the “data distribution” <italic>p</italic>
, and then changing weights so as to penalize the network for drifting away from the data distribution. In practice, we depart from <xref ref-type="disp-formula" rid="pcbi.1004554.e040">Eq 9</xref>
by using “momentum” and “weight decay” [<xref rid="pcbi.1004554.ref015" ref-type="bibr">15</xref>
], as is standard in neural-network training. Our choice of momentum and decay make this equivalent to low-pass filtering the learning signal (the right-hand sides of the equations) with an overdamped second-order system before making weight changes. Biologically, it corresponds to changes in synaptic strength having their own intrinsic dynamics.</p>
<p>Training took place in “epochs.” Data in each epoch consisted of 40,000 vectors: 40 trajectories of 1000 time steps apiece, each vector consisting of the current sensory response (proprioceptive and efference-copy) and the previous hidden-unit activities (“recurrent”; see <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3B</xref>
). On the initial time step, the recurrent units were set to all zeros and drove no weight changes. In order to accelerate convergence, and although biological implausible, weight changes were made on “minibatches” of 40 input vectors, each of which corresponded to the same time point, but from the 40 different trajectories. Fresh data (40 new trajectories) were generated every five epochs. Learning rates (<italic>ϵ</italic>
) also decayed across epochs. For the rEFH trained on an uncontrolled dynamical system, the total number of epochs was 120, and the decay was exponential: for the <italic>k</italic>
<sup>th</sup>
epoch, <inline-formula id="pcbi.1004554.e041"><alternatives><graphic xlink:href="pcbi.1004554.e041.jpg" id="pcbi.1004554.e041g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M41"><mml:mrow><mml:msubsup><mml:mi>ϵ</mml:mi>
<mml:mi>k</mml:mi>
<mml:mrow><mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msubsup>
<mml:mo>=</mml:mo>
<mml:msup><mml:mn>1.1</mml:mn>
<mml:mi>k</mml:mi>
</mml:msup>
<mml:msubsup><mml:mi>ϵ</mml:mi>
<mml:mn>0</mml:mn>
<mml:mrow><mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msubsup>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. For the case of controlled dynamics, the network was trained for 1200 epochs, with the reciprocal learning rate growing according to a sigmoidal function: <inline-formula id="pcbi.1004554.e042"><alternatives><graphic xlink:href="pcbi.1004554.e042.jpg" id="pcbi.1004554.e042g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M42"><mml:mrow><mml:msubsup><mml:mi>ϵ</mml:mi>
<mml:mi>k</mml:mi>
<mml:mrow><mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msubsup>
<mml:mo>=</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:mrow><mml:mo>(</mml:mo>
<mml:mrow><mml:mfrac><mml:mrow><mml:mn>1000</mml:mn>
</mml:mrow>
<mml:mrow><mml:mn>1</mml:mn>
<mml:mo>+</mml:mo>
<mml:mtext>exp</mml:mtext>
<mml:mo>{</mml:mo>
<mml:mo>−</mml:mo>
<mml:mi>k</mml:mi>
<mml:mo>/</mml:mo>
<mml:mn>8</mml:mn>
<mml:mo>+</mml:mo>
<mml:mn>7.5</mml:mn>
<mml:mo>}</mml:mo>
</mml:mrow>
</mml:mfrac>
<mml:mspace width="1pt"></mml:mspace>
<mml:mo>+</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:mfrac><mml:mn>1</mml:mn>
<mml:mn>2</mml:mn>
</mml:mfrac>
</mml:mrow>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:msubsup><mml:mi>ϵ</mml:mi>
<mml:mn>0</mml:mn>
<mml:mrow><mml:mo>−</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
</mml:msubsup>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. (The numbers were chosen so that the sigmoid approximately matches the exponential growth for the first 120 epochs, although their exact values are not critical.)</p>
</sec>
<sec id="sec025"><title>Testing</title>
<p>Filtering was tested on a new set of 40 trajectories (40,000 vectors). At each time step, the current “sensory” (proprioception and efference copy) and recurrent responses were fed forward to the hidden layer of the network, as in training. Unlike training, however, no samples were taken from this vector of means; instead, the real-valued vector was itself returned as the recurrent response. This is equivalent to taking several (∼ 15) samples and averaging [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
]; the means themselves were used to simplify presentation of the results, since they correspond to the maximum achievable performance of the network.</p>
<p>Formally, the solution to the filtering problem is the optimal posterior distribution over the current stimulus location, given all the observations up to this point in time: <inline-formula id="pcbi.1004554.e043"><alternatives><graphic xlink:href="pcbi.1004554.e043.jpg" id="pcbi.1004554.e043g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M43"><mml:mrow><mml:mi>p</mml:mi>
<mml:mo stretchy="false">(</mml:mo>
<mml:msub><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi mathvariant="bold-italic">θ</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:msubsup><mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
</mml:msubsup>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo stretchy="false">)</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. For the controlled dynamical system, we also ask about the posterior distribution over the controls, <inline-formula id="pcbi.1004554.e044"><alternatives><graphic xlink:href="pcbi.1004554.e044.jpg" id="pcbi.1004554.e044g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M44"><mml:mrow><mml:mi>p</mml:mi>
<mml:mo stretchy="false">(</mml:mo>
<mml:msub><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>u</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo>|</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo stretchy="false">)</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
, since they are observed only noisily at each time step. We discuss optimality below, but note here that in our case these distributions are Gaussian, so their only non-zero cumulants are mean and covariance. Generically, proving optimality of the harmonium would require showing that both these cumulants can be recovered from its hidden units at every point in time; but in the present case it is only necessary to decode the posterior mean, since it is impossible for the network to keep track of the mean without also keeping track of the covariance: incorrect estimates of the latter would result in mis-weighting of the relative reliability of current sensory information and current filter estimate, resulting in suboptimal inference of the mean at the next time step.</p>
<p>Decoding the rEFH’s hidden units exploits a trick [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
]. The representational space of the hidden units is obscure; therefore, the hidden unit activities (a real-valued vector) are passed back down through the network, i.e. into the space of the inputs. Here, the optimal decoding scheme is known: it is the center of mass of each noisy hill of activity [<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
]. This decoder was applied to hidden units at each time step, for each of the 40 testing trajectories, from which errors from the actual joint angle and control input were computed.</p>
</sec>
</sec>
<sec id="sec026"><title>The optimal filtering distribution</title>
<p>For the graphical models in Figs <xref ref-type="fig" rid="pcbi.1004554.g001">1A</xref>
and <xref ref-type="fig" rid="pcbi.1004554.g003">3A</xref>
, the solution to the filtering problem can be assimilated to a variant on the Kalman filter, and therefore computed in closed form. This is because, although the emission <inline-formula id="pcbi.1004554.e045"><alternatives><graphic xlink:href="pcbi.1004554.e045.jpg" id="pcbi.1004554.e045g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M45"><mml:mrow><mml:mi>p</mml:mi>
<mml:mo stretchy="false">(</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msub><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>θ</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
</mml:msub>
<mml:mo stretchy="false">)</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
is not a Gaussian distribution over <inline-formula id="pcbi.1004554.e046"><alternatives><graphic xlink:href="pcbi.1004554.e046.jpg" id="pcbi.1004554.e046g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M46"><mml:msubsup><mml:mi>r</mml:mi>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:math>
</alternatives>
</inline-formula>
, it is a Gaussian function of <bold><italic>θ</italic>
</bold>
<sub><italic>t</italic>
</sub>
[<xref rid="pcbi.1004554.ref004" ref-type="bibr">4</xref>
, <xref rid="pcbi.1004554.ref007" ref-type="bibr">7</xref>
] (i.e., the likelihood is an unnormalized Gaussian over <bold><italic>θ</italic>
</bold>
<sub><italic>t</italic>
</sub>
)—or more precisely, of <italic>C</italic>
<bold><italic>θ</italic>
</bold>
<sub><italic>t</italic>
</sub>
, with <italic>C</italic>
the observation matrix (see <xref ref-type="disp-formula" rid="pcbi.1004554.e018">Eq 4</xref>
)—and this is the critical requirement for the derivation of Kalman’s recursive solution. The resulting modification is small: Where the emission variance and the (Gaussian-distributed) emission appear in the standard KF equations, we substitute, respectively, the scaled tuning-curve width, <inline-formula id="pcbi.1004554.e047"><alternatives><graphic xlink:href="pcbi.1004554.e047.jpg" id="pcbi.1004554.e047g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M47"><mml:mrow><mml:msubsup><mml:mi>σ</mml:mi>
<mml:mrow><mml:mtext>tc</mml:mtext>
</mml:mrow>
<mml:mn>2</mml:mn>
</mml:msubsup>
<mml:mo>/</mml:mo>
<mml:mstyle displaystyle="true"><mml:msub><mml:mo>∑</mml:mo>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mrow><mml:mspace width="1pt"></mml:mspace>
<mml:msubsup><mml:mi>r</mml:mi>
<mml:mi>i</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:mrow>
</mml:mstyle>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
, and the center of mass of the population response, <inline-formula id="pcbi.1004554.e048"><alternatives><graphic xlink:href="pcbi.1004554.e048.jpg" id="pcbi.1004554.e048g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M48"><mml:mrow><mml:mstyle displaystyle="true"><mml:msub><mml:mo>∑</mml:mo>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mrow><mml:msub><mml:mi>ξ</mml:mi>
<mml:mi>i</mml:mi>
</mml:msub>
</mml:mrow>
</mml:mstyle>
<mml:msubsup><mml:mi>r</mml:mi>
<mml:mi>i</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>/</mml:mo>
<mml:mstyle displaystyle="true"><mml:msub><mml:mo>∑</mml:mo>
<mml:mi>i</mml:mi>
</mml:msub>
<mml:mrow><mml:msubsup><mml:mi>r</mml:mi>
<mml:mi>i</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:mrow>
</mml:mstyle>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
[<xref rid="pcbi.1004554.ref016" ref-type="bibr">16</xref>
].</p>
<p>The same applies, <italic>mutatis mutandis</italic>
, to the controls. In fact, the “controlled” case provides no additional complexity, since it corresponds to an uncontrolled third-order system (since the control has its own dynamics) whose state <bold>X</bold>
<sub><italic>t</italic>
</sub>
is the concatenation of <bold>ϴ</bold>
<sub><italic>t</italic>
</sub>
and <italic>U</italic>
<sub>t</sub>
:
<disp-formula id="pcbi.1004554.e049"><alternatives><graphic xlink:href="pcbi.1004554.e049.jpg" id="pcbi.1004554.e049g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M49"><mml:mrow><mml:mi>p</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">x</mml:mi>
</mml:mrow>
<mml:mrow><mml:mi>t</mml:mi>
<mml:mo>+</mml:mo>
<mml:mn>1</mml:mn>
</mml:mrow>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>|</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">x</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>=</mml:mo>
<mml:mi mathvariant="script">N</mml:mi>
<mml:mrow><mml:mo>(</mml:mo>
<mml:mo>Γ</mml:mo>
<mml:msubsup><mml:mrow><mml:mi mathvariant="bold">x</mml:mi>
</mml:mrow>
<mml:mi>t</mml:mi>
<mml:mrow></mml:mrow>
</mml:msubsup>
<mml:mo>+</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">μ</mml:mi>
<mml:mi>x</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msub><mml:mo>Σ</mml:mo>
<mml:mi>x</mml:mi>
</mml:msub>
<mml:mo>)</mml:mo>
</mml:mrow>
<mml:mo>,</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
<label>(10)</label>
</disp-formula>
with
<disp-formula id="pcbi.1004554.e050"><alternatives><graphic xlink:href="pcbi.1004554.e050.jpg" id="pcbi.1004554.e050g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M50"><mml:mrow><mml:mo>Γ</mml:mo>
<mml:mo>:</mml:mo>
<mml:mo>=</mml:mo>
<mml:mo>[</mml:mo>
<mml:mtable><mml:mtr><mml:mtd><mml:mn>1</mml:mn>
</mml:mtd>
<mml:mtd><mml:mo>Δ</mml:mo>
</mml:mtd>
<mml:mtd><mml:mn>0</mml:mn>
</mml:mtd>
</mml:mtr>
<mml:mtr><mml:mtd><mml:mrow><mml:mo>−</mml:mo>
<mml:mfrac><mml:mi>k</mml:mi>
<mml:mi>m</mml:mi>
</mml:mfrac>
<mml:mo>Δ</mml:mo>
</mml:mrow>
</mml:mtd>
<mml:mtd><mml:mrow><mml:mn>1</mml:mn>
<mml:mo>−</mml:mo>
<mml:mfrac><mml:mi>c</mml:mi>
<mml:mi>m</mml:mi>
</mml:mfrac>
<mml:mo>Δ</mml:mo>
</mml:mrow>
</mml:mtd>
<mml:mtd><mml:mfrac><mml:mo>Δ</mml:mo>
<mml:mi>m</mml:mi>
</mml:mfrac>
</mml:mtd>
</mml:mtr>
<mml:mtr><mml:mtd><mml:mn>0</mml:mn>
</mml:mtd>
<mml:mtd><mml:mn>0</mml:mn>
</mml:mtd>
<mml:mtd><mml:mi>α</mml:mi>
</mml:mtd>
</mml:mtr>
</mml:mtable>
<mml:mo>]</mml:mo>
<mml:mo>,</mml:mo>
<mml:mspace width="18.06749pt"></mml:mspace>
<mml:msub><mml:mi mathvariant="bold-italic">μ</mml:mi>
<mml:mi>x</mml:mi>
</mml:msub>
<mml:mo>:</mml:mo>
<mml:mo>=</mml:mo>
<mml:mo>[</mml:mo>
<mml:mtable><mml:mtr><mml:mtd><mml:msub><mml:mi mathvariant="bold-italic">μ</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msub>
</mml:mtd>
</mml:mtr>
<mml:mtr><mml:mtd><mml:msub><mml:mi>μ</mml:mi>
<mml:mi>u</mml:mi>
</mml:msub>
</mml:mtd>
</mml:mtr>
</mml:mtable>
<mml:mo>]</mml:mo>
<mml:mo>,</mml:mo>
<mml:mspace width="18.06749pt"></mml:mspace>
<mml:msub><mml:mo>Σ</mml:mo>
<mml:mi>x</mml:mi>
</mml:msub>
<mml:mo>:</mml:mo>
<mml:mo>=</mml:mo>
<mml:mo>[</mml:mo>
<mml:mtable><mml:mtr><mml:mtd><mml:msub><mml:mo>Σ</mml:mo>
<mml:mi>θ</mml:mi>
</mml:msub>
</mml:mtd>
<mml:mtd><mml:mn>0</mml:mn>
</mml:mtd>
</mml:mtr>
<mml:mtr><mml:mtd><mml:mn>0</mml:mn>
</mml:mtd>
<mml:mtd><mml:msubsup><mml:mi>σ</mml:mi>
<mml:mi>u</mml:mi>
<mml:mn>2</mml:mn>
</mml:msubsup>
</mml:mtd>
</mml:mtr>
</mml:mtable>
<mml:mo>]</mml:mo>
<mml:mo>.</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</disp-formula>
In both cases, then, the posterior (filtering) distribution over the state is always Gaussian, so at every time step, one computes the posterior mean and covariance, which can be expressed in terms of the filtering distribution at the previous time step, and of the current sensory information. A full derivation appears in <xref ref-type="supplementary-material" rid="pcbi.1004554.s001">S1 Text</xref>
.</p>
<p><xref ref-type="disp-formula" rid="pcbi.1004554.e049">Eq 10</xref>
ignores some independence statements asserted by the graph of <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3A</xref>
. In fact, an EM algorithm that accounts for them can be derived; but in our experiments, this algorithm does not achieve superior results to the “agnostic” version that tries to learn unconstrained versions of Γ, <bold><italic>μ</italic>
</bold>
<sub><italic>x</italic>
</sub>
, and Σ<sub><italic>x</italic>
</sub>
. Therefore, results for EM<sup>3</sup>
throughout use the unconstrained version of the algorithm.</p>
<sec id="sec027"><title>Benchmark models</title>
<p>Error statistics for the rEFH are compared to those from four types of model. It is simplest to think of all four types using the same filtering algorithm—the KF, modified as described to account for the Poisson emissions—but running that filter on different generative models for the observed data (<inline-formula id="pcbi.1004554.e051"><alternatives><graphic xlink:href="pcbi.1004554.e051.jpg" id="pcbi.1004554.e051g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M51"><mml:mrow><mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>r</mml:mi>
</mml:mstyle>
<mml:mi>t</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
).
<list list-type="bullet"><list-item><p><bold>PROP and EfCp</bold>
: In the simplest benchmark model, joint-angle (PROP) and control (EfCp) estimates are made simply via the center of mass on the current “sensory” population. This is the optimal decoder for populations of smoothly tiled, Gaussian-tuned, Poisson neurons [<xref rid="pcbi.1004554.ref017" ref-type="bibr">17</xref>
], under the assumption of independence through time (no dynamics). It can also be thought of as a Kalman filter applied to a generative model with infinite state-transition covariance, Σ<sub><italic>x</italic>
</sub>
. (This severs the horizontal connections in Figs <xref ref-type="fig" rid="pcbi.1004554.g001">1A</xref>
and <xref ref-type="fig" rid="pcbi.1004554.g003">3A</xref>
.)</p>
</list-item>
<list-item><p><bold>OPT</bold>
: The “optimal” model runs the KF on the true generative model for the data, i.e., using the true parameters, <inline-formula id="pcbi.1004554.e052"><alternatives><graphic xlink:href="pcbi.1004554.e052.jpg" id="pcbi.1004554.e052g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M52"><mml:mrow><mml:mo>{</mml:mo>
<mml:mi>A</mml:mi>
<mml:mo>,</mml:mo>
<mml:mo>b</mml:mo>
<mml:mo>,</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:mi>C</mml:mi>
<mml:mo>,</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:mi>α</mml:mi>
<mml:mo>,</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:mi>h</mml:mi>
<mml:mo>,</mml:mo>
<mml:msub><mml:mo>∑</mml:mo>
<mml:mi>θ</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mi>σ</mml:mi>
<mml:mi>u</mml:mi>
<mml:mn>2</mml:mn>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">ν</mml:mi>
<mml:mn>0</mml:mn>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:mtext></mml:mtext>
<mml:msub><mml:mi mathvariant="bold-italic">ϒ</mml:mi>
<mml:mn>0</mml:mn>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
.</p>
</list-item>
<list-item><p><bold>EM<sup><italic>n</italic>
</sup>
</bold>
: The rEFH was not, of course, given its parameters, but had to learn them, and only from the noisy population responses, <inline-formula id="pcbi.1004554.e053"><alternatives><graphic xlink:href="pcbi.1004554.e053.jpg" id="pcbi.1004554.e053g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M53"><mml:mrow><mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mi>T</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
, <inline-formula id="pcbi.1004554.e054"><alternatives><graphic xlink:href="pcbi.1004554.e054.jpg" id="pcbi.1004554.e054g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M54"><mml:mrow><mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>u</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mi>T</mml:mi>
<mml:mi>u</mml:mi>
</mml:msubsup>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
. A useful point of comparison, then, is the performance of a (Kalman) filter that has learned those parameters from the same noisy data, but following a learning procedure that is known to be optimal. For our linear, time-invariant systems, such learning rules can be derived: they are an implementation of expectation-maximization (EM), an algorithm that guarantees convergence to at least local optima. Applying EM to linear dynamical systems requires both a forward pass through the data, filtering, <italic>and</italic>
a backward pass for smoothing, i.e., computing the probability of the hidden state given the observations <italic>for all time</italic>
. These equations are derived in <xref ref-type="supplementary-material" rid="pcbi.1004554.s001">S1 Text</xref>
. Likewise, the algorithm must be told the order (number of states) of the latent dynamical system. Which order it was told is indicated by a superscript (e.g., EM<sup>2</sup>
). We emphasize that access to the backward pass of observations, and knowledge of the order of the latent dynamics, are advantages the EM-trained models (EM<sup><italic>n</italic>
</sup>
) enjoyed over the harmonium.</p>
</list-item>
<list-item><p><bold>OBS</bold>
: For the controlled system (Controlled dynamical system), one would also like to know how useful knowledge of the controls is to inference of the joint angle. The optimal model that ignores controls, again under the assumption of linear dynamics, can be constructed by training on fully observed data, where the model parameters <inline-formula id="pcbi.1004554.e055"><alternatives><graphic xlink:href="pcbi.1004554.e055.jpg" id="pcbi.1004554.e055g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M55"><mml:mrow><mml:mo>{</mml:mo>
<mml:mi>A</mml:mi>
<mml:mo>,</mml:mo>
<mml:mspace width="1pt"></mml:mspace>
<mml:mi>C</mml:mi>
<mml:mo>,</mml:mo>
<mml:msub><mml:mo>∑</mml:mo>
<mml:mi>θ</mml:mi>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:msub><mml:mi mathvariant="bold-italic">ν</mml:mi>
<mml:mn>0</mml:mn>
</mml:msub>
<mml:mo>,</mml:mo>
<mml:mtext></mml:mtext>
<mml:msub><mml:mi mathvariant="bold-italic">ϒ</mml:mi>
<mml:mn>0</mml:mn>
</mml:msub>
<mml:mo>}</mml:mo>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
are learned from <bold>ϴ</bold>
<sub>0</sub>
, …, <bold>ϴ</bold>
<sub><italic>T</italic>
</sub>
, <inline-formula id="pcbi.1004554.e056"><alternatives><graphic xlink:href="pcbi.1004554.e056.jpg" id="pcbi.1004554.e056g" mimetype="image" position="anchor" orientation="portrait"></graphic>
<mml:math id="M56"><mml:mrow><mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mn>0</mml:mn>
<mml:mi>θ</mml:mi>
</mml:msubsup>
<mml:mo>,</mml:mo>
<mml:mo>…</mml:mo>
<mml:mo>,</mml:mo>
<mml:msubsup><mml:mstyle mathvariant="bold" mathsize="normal"><mml:mi>R</mml:mi>
</mml:mstyle>
<mml:mi>T</mml:mi>
<mml:mi>θ</mml:mi>
</mml:msubsup>
</mml:mrow>
</mml:math>
</alternatives>
</inline-formula>
, essentially via linear regression (see <xref ref-type="supplementary-material" rid="pcbi.1004554.s001">S1 Text</xref>
). This corresponds to using the generative model of <xref ref-type="fig" rid="pcbi.1004554.g001">Fig 1A</xref>
, even though the true data were generated by the model of <xref ref-type="fig" rid="pcbi.1004554.g003">Fig 3A</xref>
. OBS is therefore suboptimal, but tells us how much suboptimality accrues by ignoring the controls. (We fit the parameters of OBS, rather than providing them, because the fit model can actually outperform one based on the true dynamical-system parameters. This is because OBS can compensate to some extent for the missing controls by overestimating the state-transition noise. If the model were forced to use the true state-transition noise, but still assume zero control input, it would be worse at explaining state transitions.)</p>
</list-item>
</list>
The benchmark models also enjoyed another advantage over the rEFH. The sufficient statistics of the emission for the state are the population center of mass and the scaled tuning-curve width (or simply the scaling factor), as alluded to above. These were given to the benchmark models, rather than learned. (In the standard model, linear-Gaussian emissions with fixed variance, learning that variance with EM is straightforward [<xref rid="pcbi.1004554.ref018" ref-type="bibr">18</xref>
]. For our more complicated emission model, it is not, which is why we decided simply to provide it to the benchmark models.) The harmonium, on the other hand, had to learn what to do with the vectors of raw spike counts, <bold>r</bold>
<sup><italic>θ</italic>
</sup>
, <bold>r</bold>
<sup><italic>u</italic>
</sup>
.</p>
<p>Different training runs, like different testing sets, will yield slightly different models. Thus, for each type of model to be trained, including the rEFH, we selected the best of 20 different networks, each trained from scratch with random initialization. Then we repeated this procedure itself twelve times, and from these twelve tokens of each model type, generated the error bars for the MSEs in Figs <xref ref-type="fig" rid="pcbi.1004554.g001">1E</xref>
, <xref ref-type="fig" rid="pcbi.1004554.g003">3D</xref>
and <xref ref-type="fig" rid="pcbi.1004554.g003">3F</xref>
. Each of the twelve tokens of a particular model type was tested on a different testing set, but the <italic>same</italic>
testing set was used for matching tokens of different types (so, e.g., the fourth rEFH token was tested on the same data as the fourth EM<sup>2</sup>
token).</p>
</sec>
</sec>
<sec id="sec028"><title>Tuning analysis</title>
<p>In the section <bold>Learned receptive fields and connectivity</bold>
, in order to determine how the network has learned to solve the filtering problem, we sort hidden units by their “preferred” lags and “preferred” angles. These were computed as follows. First, we generated a new set of 40 trajectories of 1000 time steps apiece. Then we computed hidden-unit mean activities, i.e., their probability of firing (these are the same quantity because the hidden units are conditionally Bernoulli random variables). Angular positions for all 40000 time points were then discretized into 30 bins of uniform width spanning the feasible joint space.</p>
<p>For each hidden unit, the following calculation was then performed. First, the empirical mutual information (MI) was computed, according to the standard formula [<xref rid="pcbi.1004554.ref019" ref-type="bibr">19</xref>
], between the two discrete random variables: the discretized position (30 bins) and the binary (spike/no spike) hidden-unit response. Next, to reject spurious MI (which will anyway be rare, given the number of data), for each of 20 reshuffles, the unit’s response was shuffled in time and the MIs recalculated. If the unit’s unshuffled MI fell below the 95th percentile of its shuffled MIs, the unshuffled MI was set to zero. The entire procedure was then repeated with the response time-shifted forward by one step, for each of 40 steps. Finally, the “preferred” lag was selected to be the time shift for which MI was maximized. These were used to sort the receptive fields in Figs <xref ref-type="fig" rid="pcbi.1004554.g005">5A and 5B</xref>
, <xref ref-type="fig" rid="pcbi.1004554.g007">7</xref>
and <xref ref-type="fig" rid="pcbi.1004554.g008">8B</xref>
.</p>
<p>For each unit, a “lagged” tuning curve can be constructed by considering its mean responses to past (discretized) stimuli; in particular, to stimuli at that unit’s <italic>preferred</italic>
lag. These are the curves plotted as a heat map in <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5C</xref>
, where they have been sorted by the locations of the tuning curves’ peaks. The same locations were used to sort the weight matrix in <xref ref-type="fig" rid="pcbi.1004554.g008">Fig 8A</xref>
. Inverting the process, one can ask how well these tuning curves explain the receptive fields in the space of non-delayed position and velocity (<xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5A</xref>
): apply each tuning curve to each of the 40000 stimuli, delay the responses by the units’ preferred lags, and then compute receptive fields with these responses. This is how <xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5B</xref>
was constructed.</p>
<p>Finally, comparing the distribution of preferred lags (<xref ref-type="fig" rid="pcbi.1004554.g005">Fig 5D</xref>
) to the autocorrelation of the stimulus required computing the autocorrelation of a circular variable (angle). We used the angular-angular correlation measure given by Zar [<xref rid="pcbi.1004554.ref020" ref-type="bibr">20</xref>
].</p>
</sec>
</sec>
<sec sec-type="supplementary-material" id="sec029"><title>Supporting Information</title>
<supplementary-material content-type="local-data" id="pcbi.1004554.s001"><label>S1 Text</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s001.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
<supplementary-material content-type="local-data" id="pcbi.1004554.s002"><label>S2 Text</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s002.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
<supplementary-material content-type="local-data" id="pcbi.1004554.s003"><label>S3 Text</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s003.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
<supplementary-material content-type="local-data" id="pcbi.1004554.s004"><label>S1 Fig</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s004.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
<supplementary-material content-type="local-data" id="pcbi.1004554.s005"><label>S2 Fig</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s005.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
<supplementary-material content-type="local-data" id="pcbi.1004554.s006"><label>S4 Text</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s006.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
<supplementary-material content-type="local-data" id="pcbi.1004554.s007"><label>S3 Fig</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s007.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
<supplementary-material content-type="local-data" id="pcbi.1004554.s008"><label>S1 Bibliography</label>
<caption><p>(PDF)</p>
</caption>
<media xlink:href="pcbi.1004554.s008.pdf"><caption><p>Click here for additional data file.</p>
</caption>
</media>
</supplementary-material>
</sec>
</body>
<back><ack><p>Some of the EFHs were trained using Tesla K40 GPUs, the generous donation of the Nvidia Corporation.</p>
</ack>
<ref-list><title>References</title>
<ref id="pcbi.1004554.ref001"><label>1</label>
<mixed-citation publication-type="book"><name><surname>Földiák</surname>
<given-names>P</given-names>
</name>
. <chapter-title>The “Ideal Homunculus”: Statistical Inference from Neural Population Responses</chapter-title>
In: <name><surname>Eeckman</surname>
<given-names>FH</given-names>
</name>
, <name><surname>Bower</surname>
<given-names>JM</given-names>
</name>
, editors. <source>Computation and Neural Systems</source>
. Norwell, MA: <publisher-loc>Norwell, MA</publisher-loc>
: <publisher-name>Kluwer Academic Publishers</publisher-name>
; <year>1993</year>
p. <fpage>55</fpage>
–<lpage>60</lpage>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref002"><label>2</label>
<mixed-citation publication-type="journal"><name><surname>Ma</surname>
<given-names>WJ</given-names>
</name>
, <name><surname>Beck</surname>
<given-names>JM</given-names>
</name>
, <name><surname>Latham</surname>
<given-names>PE</given-names>
</name>
, <name><surname>Pouget</surname>
<given-names>A</given-names>
</name>
. <article-title>Bayesian Inference with Probabilistic Population Codes</article-title>
. <source>Nature Neuroscience</source>
. <year>2006</year>
;<volume>9</volume>
:<fpage>1423</fpage>
–<lpage>1438</lpage>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref003"><label>3</label>
<mixed-citation publication-type="other">Welling M, Rosen-Zvi M, Hinton GE. Exponential Family Harmoniums with an Application to Information Retrieval. In: Advances in Neural Information Processing Systems 17: Proceedings of the 2004 Conference; 2005. p. 1481–1488.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref004"><label>4</label>
<mixed-citation publication-type="journal"><name><surname>Makin</surname>
<given-names>JG</given-names>
</name>
, <name><surname>Fellows</surname>
<given-names>MR</given-names>
</name>
, <name><surname>Sabes</surname>
<given-names>PN</given-names>
</name>
. <article-title>Learning Multisensory Integration and Coordinate Transformation via Density Estimation</article-title>
. <source>PLoS Computational Biology</source>
. <year>2013</year>
;<volume>9</volume>
(<issue>4</issue>
):<fpage>1</fpage>
–<lpage>17</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1371/journal.pcbi.1003035">10.1371/journal.pcbi.1003035</ext-link>
</comment>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref005"><label>5</label>
<mixed-citation publication-type="journal"><name><surname>van Beers</surname>
<given-names>RJ</given-names>
</name>
, <name><surname>Sittig</surname>
<given-names>A</given-names>
</name>
, <name><surname>van Der Gon</surname>
<given-names>JJD</given-names>
</name>
. <article-title>Integration of Proprioceptive and Visual Position-Information: An Experimentally Supported Model</article-title>
. <source>Journal of Neurophysiology</source>
. <year>1999</year>
;<volume>81</volume>
:<fpage>1355</fpage>
–<lpage>1364</lpage>
. <pub-id pub-id-type="pmid">10085361</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref006"><label>6</label>
<mixed-citation publication-type="journal"><name><surname>Ernst</surname>
<given-names>MO</given-names>
</name>
, <name><surname>Banks</surname>
<given-names>MS</given-names>
</name>
. <article-title>Humans Integrate Visual and Haptic Information in a Statistically Optimal Fashion</article-title>
. <source>Nature</source>
. <year>2002</year>
;<volume>415</volume>
(<issue>January</issue>
):<fpage>429</fpage>
–<lpage>433</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1038/415429a">10.1038/415429a</ext-link>
</comment>
<pub-id pub-id-type="pmid">11807554</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref007"><label>7</label>
<mixed-citation publication-type="journal"><name><surname>Beck</surname>
<given-names>JM</given-names>
</name>
, <name><surname>Latham</surname>
<given-names>PE</given-names>
</name>
, <name><surname>Pouget</surname>
<given-names>A</given-names>
</name>
. <article-title>Marginalization in Neural Circuits with Divisive Normalization</article-title>
. <source>Journal of Neuroscience</source>
. <year>2011</year>
<month>10</month>
;<volume>31</volume>
(<issue>43</issue>
):<fpage>15310</fpage>
–<lpage>9</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1523/JNEUROSCI.1706-11.2011">10.1523/JNEUROSCI.1706-11.2011</ext-link>
</comment>
<pub-id pub-id-type="pmid">22031877</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref008"><label>8</label>
<mixed-citation publication-type="journal"><name><surname>Boerlin</surname>
<given-names>M</given-names>
</name>
, <name><surname>Denève</surname>
<given-names>S</given-names>
</name>
. <article-title>Spike-Based Population Coding and Working Memory</article-title>
. <source>PLoS Computational Biology</source>
. <year>2011</year>
<month>2</month>
;<volume>7</volume>
(<issue>2</issue>
):<fpage>e1001080</fpage>
<comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1371/journal.pcbi.1001080">10.1371/journal.pcbi.1001080</ext-link>
</comment>
<pub-id pub-id-type="pmid">21379319</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref009"><label>9</label>
<mixed-citation publication-type="other">Makin JG, Sabes PN. Sensory Integration and Density Estimation. Advances in Neural Information Processing Systems 27: Proceedings of the 2014 Conference. 2015;p. 1–9.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref010"><label>10</label>
<mixed-citation publication-type="journal"><name><surname>Wolpert</surname>
<given-names>DM</given-names>
</name>
, <name><surname>Goodbody</surname>
<given-names>SJ</given-names>
</name>
, <name><surname>Husain</surname>
<given-names>M</given-names>
</name>
. <article-title>Maintaining Internal Representations: the Role of the Human Superior Parietal Lobe</article-title>
. <source>Nature Neuroscience</source>
. <year>1998</year>
;<volume>1</volume>
:<fpage>529</fpage>
–<lpage>533</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1038/2245">10.1038/2245</ext-link>
</comment>
<pub-id pub-id-type="pmid">10196553</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref011"><label>11</label>
<mixed-citation publication-type="journal"><name><surname>Mulliken</surname>
<given-names>GH</given-names>
</name>
, <name><surname>Musallam</surname>
<given-names>S</given-names>
</name>
, <name><surname>Andersen</surname>
<given-names>RA</given-names>
</name>
. <article-title>Decoding trajectories from posterior parietal cortex ensembles</article-title>
. <source>Journal of Neuroscience</source>
. <year>2008</year>
<month>11</month>
;<volume>28</volume>
(<issue>48</issue>
):<fpage>12913</fpage>
–<lpage>26</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1523/JNEUROSCI.1463-08.2008">10.1523/JNEUROSCI.1463-08.2008</ext-link>
</comment>
<pub-id pub-id-type="pmid">19036985</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref012"><label>12</label>
<mixed-citation publication-type="journal"><name><surname>Hinton</surname>
<given-names>GE</given-names>
</name>
, <name><surname>Osindero</surname>
<given-names>S</given-names>
</name>
, <name><surname>Teh</surname>
<given-names>YW</given-names>
</name>
. <article-title>A Fast Learning Algorithm for Deep Belief Nets</article-title>
. <source>Neural Computation</source>
. <year>2006</year>
;<volume>18</volume>
:<fpage>1527</fpage>
–<lpage>1554</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1162/neco.2006.18.7.1527">10.1162/neco.2006.18.7.1527</ext-link>
</comment>
<pub-id pub-id-type="pmid">16764513</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref013"><label>13</label>
<mixed-citation publication-type="journal"><name><surname>Hinton</surname>
<given-names>GE</given-names>
</name>
. <article-title>Training Products of Experts by Minimizing Contrastive Divergence</article-title>
. <source>Neural Computation</source>
. <year>2002</year>
;<volume>14</volume>
:<fpage>1771</fpage>
–<lpage>1800</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1162/089976602760128018">10.1162/089976602760128018</ext-link>
</comment>
<pub-id pub-id-type="pmid">12180402</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref014"><label>14</label>
<mixed-citation publication-type="book"><name><surname>Scherpen</surname>
<given-names>JMA</given-names>
</name>
. <chapter-title>Chapter 4: Balanced Realizations, Model Order Reduction, and the Hankel Operator</chapter-title>
In: <name><surname>Levine</surname>
<given-names>WS</given-names>
</name>
, editor. <source>The Control Handbook</source>
. <edition>2nd ed</edition>
<publisher-name>CRC Press</publisher-name>
; <year>2015</year>
p. <fpage>4–1</fpage>
–<lpage>4–24</lpage>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref015"><label>15</label>
<mixed-citation publication-type="book"><name><surname>Hinton</surname>
<given-names>GE</given-names>
</name>
. <source>A Practical Guide to Training Restricted Boltzmann Machines</source>
. <publisher-loc>Toronto</publisher-loc>
: <publisher-name>University of Toronto</publisher-name>
; <year>2010</year>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref016"><label>16</label>
<mixed-citation publication-type="journal"><name><surname>Makin</surname>
<given-names>JG</given-names>
</name>
, <name><surname>Fellows</surname>
<given-names>MR</given-names>
</name>
, <name><surname>Sabes</surname>
<given-names>PN</given-names>
</name>
. <article-title>Learning Multisensory Integration and Coordinate Transformation via Density Estimation—Supporting Material</article-title>
. <source>PLoS Computational Biology</source>
. <year>2013</year>
;<volume>9</volume>
(<issue>4</issue>
):<fpage>1</fpage>
–<lpage>9</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1371/journal.pcbi.1003035">10.1371/journal.pcbi.1003035</ext-link>
</comment>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref017"><label>17</label>
<mixed-citation publication-type="book"><name><surname>Dayan</surname>
<given-names>P</given-names>
</name>
, <name><surname>Abbott</surname>
<given-names>LF</given-names>
</name>
. <source>Theoretical Neuroscience</source>
. <publisher-name>MIT Press</publisher-name>
; <year>2001</year>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref018"><label>18</label>
<mixed-citation publication-type="book"><name><surname>Ghahramani</surname>
<given-names>Z</given-names>
</name>
, <name><surname>Hinton</surname>
<given-names>GE</given-names>
</name>
. <source>Parameter Estimation for Linear Dynamical Systems</source>
. <publisher-name>University of Toronto</publisher-name>
; <year>1996</year>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref019"><label>19</label>
<mixed-citation publication-type="book"><name><surname>Cover</surname>
<given-names>TM</given-names>
</name>
, <name><surname>Thomas</surname>
<given-names>JA</given-names>
</name>
. <source>Elements of Information Theory</source>
. <publisher-name>Wiley</publisher-name>
; <year>2006</year>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref020"><label>20</label>
<mixed-citation publication-type="book"><name><surname>Zar</surname>
<given-names>JH</given-names>
</name>
. <source>Biostatistical Analysis</source>
. <publisher-name>Prentice Hall</publisher-name>
; <year>1999</year>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref021"><label>21</label>
<mixed-citation publication-type="journal"><name><surname>Andersen</surname>
<given-names>RA</given-names>
</name>
, <name><surname>Snyder</surname>
<given-names>LH</given-names>
</name>
, <name><surname>Bradley</surname>
<given-names>DC</given-names>
</name>
, <name><surname>Xing</surname>
<given-names>J</given-names>
</name>
. <article-title>Multimodal Representation of Space in the Posterior Parietal Cortex and its Use in Planning Movements</article-title>
. <source>Annual Review of Neuroscience</source>
. <year>1997</year>
;<volume>20</volume>
:<fpage>303</fpage>
–<lpage>330</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1146/annurev.neuro.20.1.303">10.1146/annurev.neuro.20.1.303</ext-link>
</comment>
<pub-id pub-id-type="pmid">9056716</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref022"><label>22</label>
<mixed-citation publication-type="journal"><name><surname>Egger</surname>
<given-names>SW</given-names>
</name>
, <name><surname>Britten</surname>
<given-names>KH</given-names>
</name>
. <article-title>Linking sensory neurons to visually guided behavior: relating MST activity to steering in a virtual environment</article-title>
. <source>Visual neuroscience</source>
. <year>2013</year>
<month>11</month>
;<volume>30</volume>
(<issue>5–6</issue>
):<fpage>315</fpage>
–<lpage>30</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1017/S0952523813000412">10.1017/S0952523813000412</ext-link>
</comment>
<pub-id pub-id-type="pmid">24171813</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref023"><label>23</label>
<mixed-citation publication-type="journal"><name><surname>Kuenzle</surname>
<given-names>H</given-names>
</name>
. <article-title>Cortico-Cortical Efferents of Primary Motor and Somatosensory Regions of the Cerebral Cortex in Macaca Fascicularis</article-title>
. <source>Neuroscience</source>
. <year>1978</year>
;<volume>3</volume>
:<fpage>25</fpage>
–<lpage>39</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1016/0306-4522(78)90151-3">10.1016/0306-4522(78)90151-3</ext-link>
</comment>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref024"><label>24</label>
<mixed-citation publication-type="journal"><name><surname>Ghosh</surname>
<given-names>S</given-names>
</name>
, <name><surname>Brinkman</surname>
<given-names>C</given-names>
</name>
, <name><surname>Porter</surname>
<given-names>R</given-names>
</name>
. <article-title>A Quantitative Study of the Distribution of Neurons Projecting to the Precentral Motor Cortex in the Monkey (M. Fascicularis)</article-title>
. <source>The Journal of Comparative Neurology</source>
. <year>1987</year>
;<volume>259</volume>
:<fpage>424</fpage>
–<lpage>444</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1002/cne.902590309">10.1002/cne.902590309</ext-link>
</comment>
<pub-id pub-id-type="pmid">3584565</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref025"><label>25</label>
<mixed-citation publication-type="journal"><name><surname>Burchinskaya</surname>
<given-names>LF</given-names>
</name>
. <article-title>Neuronal Composition and Interneuronal Connection of Area 5 in the Cat Parietal Association Cortex</article-title>
. <source>Neirofiziologiya</source>
. <year>1979</year>
;<volume>11</volume>
(<issue>1</issue>
):<fpage>35</fpage>
–<lpage>42</lpage>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref026"><label>26</label>
<mixed-citation publication-type="journal"><name><surname>Carracedo</surname>
<given-names>LM</given-names>
</name>
, <name><surname>Kjeldsen</surname>
<given-names>H</given-names>
</name>
, <name><surname>Cunnington</surname>
<given-names>L</given-names>
</name>
, <name><surname>Jenkins</surname>
<given-names>a</given-names>
</name>
, <name><surname>Schofield</surname>
<given-names>I</given-names>
</name>
, <name><surname>Cunningham</surname>
<given-names>MO</given-names>
</name>
, <etal>et al</etal>
<article-title>A Neocortical Delta Rhythm Facilitates Reciprocal Interlaminar Interactions via Nested Theta Rhythms</article-title>
. <source>Journal of Neuroscience</source>
. <year>2013</year>
;<volume>33</volume>
(<issue>26</issue>
):<fpage>10750</fpage>
–<lpage>10761</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1523/JNEUROSCI.0735-13.2013">10.1523/JNEUROSCI.0735-13.2013</ext-link>
</comment>
<pub-id pub-id-type="pmid">23804097</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref027"><label>27</label>
<mixed-citation publication-type="journal"><name><surname>Rao</surname>
<given-names>RPN</given-names>
</name>
, <name><surname>Ballard</surname>
<given-names>DH</given-names>
</name>
. <article-title>Dynamic Model of Visual Recognition Predicts Neural Response Properties in the Visual Cortex</article-title>
. <source>Neural Computation</source>
. <year>1997</year>
<month>5</month>
;<volume>9</volume>
(<issue>4</issue>
):<fpage>721</fpage>
–<lpage>63</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1162/neco.1997.9.4.721">10.1162/neco.1997.9.4.721</ext-link>
</comment>
<pub-id pub-id-type="pmid">9161021</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref028"><label>28</label>
<mixed-citation publication-type="journal"><name><surname>Denève</surname>
<given-names>S</given-names>
</name>
, <name><surname>Duhamel</surname>
<given-names>JR</given-names>
</name>
, <name><surname>Pouget</surname>
<given-names>A</given-names>
</name>
. <article-title>Optimal Sensorimotor Integration in Recurrent Cortical Networks: A Neural Implementation of Kalman Filters</article-title>
. <source>Journal of Neuroscience</source>
. <year>2007</year>
;<volume>27</volume>
(<issue>21</issue>
):<fpage>5744</fpage>
–<lpage>5756</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1523/JNEUROSCI.3985-06.2007">10.1523/JNEUROSCI.3985-06.2007</ext-link>
</comment>
<pub-id pub-id-type="pmid">17522318</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref029"><label>29</label>
<mixed-citation publication-type="journal"><name><surname>Huys</surname>
<given-names>QJM</given-names>
</name>
, <name><surname>Zemel</surname>
<given-names>RS</given-names>
</name>
, <name><surname>Natarajan</surname>
<given-names>R</given-names>
</name>
, <name><surname>Dayan</surname>
<given-names>P</given-names>
</name>
. <article-title>Fast Population Coding</article-title>
. <source>Neural Computation</source>
. <year>2007</year>
;<volume>19</volume>
:<fpage>404</fpage>
–<lpage>441</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1162/neco.2007.19.2.404">10.1162/neco.2007.19.2.404</ext-link>
</comment>
<pub-id pub-id-type="pmid">17206870</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref030"><label>30</label>
<mixed-citation publication-type="journal"><name><surname>Natarajan</surname>
<given-names>R</given-names>
</name>
, <name><surname>Huys</surname>
<given-names>QJM</given-names>
</name>
, <name><surname>Dayan</surname>
<given-names>P</given-names>
</name>
, <name><surname>Zemel</surname>
<given-names>RS</given-names>
</name>
. <article-title>Encoding and decoding spikes for dynamic stimuli</article-title>
. <source>Neural Computation</source>
. <year>2008</year>
;<volume>20</volume>
:<fpage>2325</fpage>
–<lpage>2360</lpage>
. <comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1162/neco.2008.01-07-436">10.1162/neco.2008.01-07-436</ext-link>
</comment>
<pub-id pub-id-type="pmid">18386986</pub-id>
</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref031"><label>31</label>
<mixed-citation publication-type="other">Hinton GE, Brown A. Spiking Boltzmann Machines. Advances in Neural Information Processing Systems 12: Proceedings of the 1999 Conference. 2000;12.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref032"><label>32</label>
<mixed-citation publication-type="book"><name><surname>Sutskever</surname>
<given-names>I</given-names>
</name>
, <name><surname>Hinton</surname>
<given-names>GE</given-names>
</name>
. <source>Learning Multilevel Distributed Representations for High-Dimensional Sequences</source>
. In: <publisher-name>AISTATS</publisher-name>
; <year>2007</year>
p. <fpage>1</fpage>
–<lpage>8</lpage>
.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref033"><label>33</label>
<mixed-citation publication-type="other">Sutskever I, Hinton GE, Taylor G. The Recurrent Temporal Restricted Boltzmann Machine. In: Advances in Neural Information Processing Systems 21: Proceedings of the 2008 Conference; 2009. p. 1–8.</mixed-citation>
</ref>
<ref id="pcbi.1004554.ref034"><label>34</label>
<mixed-citation publication-type="journal"><name><surname>Wolpert</surname>
<given-names>DM</given-names>
</name>
, <name><surname>Ghahramani</surname>
<given-names>Z</given-names>
</name>
, <name><surname>Jordan</surname>
<given-names>MI</given-names>
</name>
. <article-title>An Internal Model for Sensorimotor Integration</article-title>
. <source>Science</source>
. <year>1995</year>
;<volume>269</volume>
(<issue>5232</issue>
):<fpage>1880</fpage>
<comment>doi: <ext-link ext-link-type="uri" xlink:href="http://dx.doi.org/10.1126/science.7569931">10.1126/science.7569931</ext-link>
</comment>
<pub-id pub-id-type="pmid">7569931</pub-id>
</mixed-citation>
</ref>
</ref-list>
</back>
</pmc>
<affiliations><list><country><li>États-Unis</li>
</country>
<region><li>Californie</li>
</region>
</list>
<tree><country name="États-Unis"><region name="Californie"><name sortKey="Makin, Joseph G" sort="Makin, Joseph G" uniqKey="Makin J" first="Joseph G." last="Makin">Joseph G. Makin</name>
</region>
<name sortKey="Dichter, Benjamin K" sort="Dichter, Benjamin K" uniqKey="Dichter B" first="Benjamin K." last="Dichter">Benjamin K. Dichter</name>
<name sortKey="Dichter, Benjamin K" sort="Dichter, Benjamin K" uniqKey="Dichter B" first="Benjamin K." last="Dichter">Benjamin K. Dichter</name>
<name sortKey="Makin, Joseph G" sort="Makin, Joseph G" uniqKey="Makin J" first="Joseph G." last="Makin">Joseph G. Makin</name>
<name sortKey="Sabes, Philip N" sort="Sabes, Philip N" uniqKey="Sabes P" first="Philip N." last="Sabes">Philip N. Sabes</name>
<name sortKey="Sabes, Philip N" sort="Sabes, Philip N" uniqKey="Sabes P" first="Philip N." last="Sabes">Philip N. Sabes</name>
<name sortKey="Sabes, Philip N" sort="Sabes, Philip N" uniqKey="Sabes P" first="Philip N." last="Sabes">Philip N. Sabes</name>
</country>
</tree>
</affiliations>
</record>
Pour manipuler ce document sous Unix (Dilib)
EXPLOR_STEP=$WICRI_ROOT/Ticri/CIDE/explor/HapticV1/Data/Ncbi/Merge
HfdSelect -h $EXPLOR_STEP/biblio.hfd -nk 003D39 | SxmlIndent | more
Ou
HfdSelect -h $EXPLOR_AREA/Data/Ncbi/Merge/biblio.hfd -nk 003D39 | SxmlIndent | more
Pour mettre un lien sur cette page dans le réseau Wicri
{{Explor lien |wiki= Ticri/CIDE |area= HapticV1 |flux= Ncbi |étape= Merge |type= RBID |clé= PMC:4634970 |texte= Learning to Estimate Dynamical State with Probabilistic Population Codes }}
Pour générer des pages wiki
HfdIndexSelect -h $EXPLOR_AREA/Data/Ncbi/Merge/RBID.i -Sk "pubmed:26540152" \ | HfdSelect -Kh $EXPLOR_AREA/Data/Ncbi/Merge/biblio.hfd \ | NlmPubMed2Wicri -a HapticV1
![]() | This area was generated with Dilib version V0.6.23. | ![]() |