FY2016 Annual Report

Computational Neuroscience Unit

Professor Erik De Schutter

 

Abstract

We use computational, data-driven methods to study how neurons and microcircuits in the brain operate. We are interested  in how fundamental properties, such as a neuron’s morphology and its excitability, interact with one another during common neural functions like information processing or learning. Most of our models concern the cerebellum as this brain structure has a relatively simple anatomy and the physiology of its main neurons has been studied extensively, allowing for detailed modeling at many different levels of complexity.

 

1. Staff

Molecular modeling

• Iain Hepburn, Technical Staff
• Andrew Gallimore, Researcher 
• Criseida Zamora, Researcher 

Cellular modeling

• Sungho Hong, Group Leader
• Yunliang Zang, Researcher

Network modeling

• Sergio Verduzco, Researcher 
• Peter Bratby, Researcher (from November 2016)
• Mykola Medvidov, Staff Scientist (from March 2017)

Software development

• Weiliang Chen, Researcher
• Guido Klingbeil, Researcher

Visiting Researcher

• Marylka Yoe Uusisaari (till December 2016)
• Akira Takashima (from June 2016)
• Yoshiyuki Kubota (from October 2016)
• ​Jihwan Myung (from November 2016)

Research Interns

• Hanzhi Chen, Research Intern (June - August 2016)
• Yoshino Hasegawa, Research Intern (August - September 2016)
• Claus Lang, Research Intern (October - December 2016)
• Chie Taguchi, Research Intern (March 2016)

Rotation Students

• Mizuki Kato, Student (Term 1)
• Nadine Wilkuttis, Student (Term  1)
• Dongqi Han, Student (Term 2)
• Ayaka Usui, Student (Term 3)

Research Unit Administrator

• Sachie Matsuoka

​

2. Collaborations

• Theme: Cerebellar physiology, multiple themes
  • Type of collaboration: Scientific collaboration and graduate program
  • Researchers:
    • Professor M. Giugliano, University of Antwerp, Belgium
    • Professor D. Snyders, University of Antwerp, Belgium
    • Joao Couto, University of Antwerp, Belgium
    • Q. Robberecht, University of Antwerp, Belgium
• Theme: Spiking activity of monkey cerebellar neurons
  • Type of collaboration: Scientific collaboration
  • Researchers:
    • Professor H.P. Thier, University of Tübingen, Germany
    • A. Ignashchenkova, University of Tübingen, Germany
    • Dr. M. Junker, University of Tübingen, Germany
    • A. Schmigdlin, University of Tübingen, Germany
       
• Theme: Human Brain Project: simulator development
  • Type of collaboration: Scientific collaboration
  • Researchers:
    • Prof. F. Schürmann, École Polytechnique Fédérale de Lausanne, Switzerland
    • F. Delalondre, École Polytechnique Fédérale de Lausanne, Switzerland

• Theme: Molecular identification of cerebellar signaling pathways and cerebellar optogenetics
  • Type of collaboration: Scientific collaboration 
  • Researchers:
    • Professor K. Tanaka, Korea Institute for Science and Technology (KIST), Korea
       
• Theme: Modeling of effects of ethanol on the cerebellum
  • Type of collaboration: Joint research
  • Researchers:
    • Professor C.F. Valenzuela, University of New Mexico, United States of America
       
• Theme: In vivo cerebellar activity
  • Type of collaboration: Scientific collaboration
  • Researchers:
    • Professor C. De Zeeuw, Erasmus Medical Center, Rotterdam, The Netherlands
    • Professor L.W.J. Bosman, Erasmus Medical Center, Rotterdam, The Netherlands
    • Dr. M. Negrello, Erasmus Medical Center, Rotterdam, The Netherlands

• Theme: Purkinje cell morphology and physiology, modeling
  • Type of collaboration: Scientific collaboration
  • Researchers:
    • Professor M. Häusser, University College London, United Kingdom
    • Professor A. Roth, University College London, United Kingdom
    • Dr. S. Dieudonné, Ecole Normale Supérieure, Paris, France

3. Activities and Findings

3.1 STEPS software development

Accurate reaction-diusion operator splitting on tetrahedral meshes 

The Inhomogeneous SSA (ISSA) used in STEPS is perfectly accurate but also quite slow, especially for large systems. This became a real bottleneck for our research, for example simulations of stochastic calcium spikes  took more than 3 weeks to compute and this was for only a section representing about 15% of the dendritic tree (Anwar et al, 2013).  A solution would be to parallelize the software so that it can be run on multiple cores at the same time but we had to develop a method to do this efficiently. The standard ISSA cannot be parallelized effectively because diffusion events tend to dominate the simulation and this causes a large amount of  communication to handle diffusion events between tetrahedrons running on different cores.

The solution to this problem is an Operator Split approach, where the chemical reactions run according to the standard SSA but diffusion events are executed only at specific time steps, allowing for buffered communication (Figure 1). We systematically investigated the existing literature on stochastic reaction-diffusion operator splitting, but most of this applied to regular meshes and some approaches suffered from poor accuracy (Hepburn, Chen & De Schutter, 2016). We therefore developed a new operator splitting implementation for irregular meshes that enhances accuracy with minimal performance cost (Hepburn, Chen & De Schutter, 2016). The method requires a global value τ that can be computed before the simulation starts and allows for multiple reaction and diffusion events in a single time step. 

Initial parallel implementation in STEPS

Using this Operator Split approac  we developed an initial MPI-based, parallel implementation of STEPS without E-Field (Chen and De Schutter 2017) as shown in Figure 2.

Figure 2: Operator Split method leads to more effective communication of diffusion events between two processes p1 and p2 in a parallel solution. Note that the reaction SSA continues in process 2 while it is still communicating (from Chen and De Schutter, 2017).

Simulation results show that this initial implementation is capable of achieving super-linear speedup for simulations with reasonable molecule density and mesh quality and a balanced loading; the super-linearity is due to more effective use of processor memory caching. In the best scenario, a parallel simulation with 2,000 processes runs more than 3,600 times faster than its serial SSA counterpart, and achieves more than 20-fold speedup relative to parallel simulation with 100 processes (Figure 3). In a more realistic scenario with dynamic calcium influx and data recording, the parallel simulation with 1,000 processes and no load balancing is still 500 times faster than the conventional serial SSA simulation. 

3.2 Molecular mechanisms of synaptic plasticity

Improved model of stochastic cerebellar LTD induction
Models should evolve as new experimental data becomes available. When we developed the Antunes and De Schutter (2012) model it was unclear how PKC contributed to activation of Raf, the first step in activation of the ERK-based positive feedback loop. Therefore a placeholder molecule ‘Raf-activator’ was used in this model. Subsequently, Yamamoto et al. (2012) showed that the actual mechanism occurs through Raf Kinase Inhibitory Protein (RKIP). RKIP binds to Raf and inactivates it, phosphorylation of RKIP by PKC causes the complex to dissociate and releases active Raf. To model this interaction we borrowed from a previous model of RKIP activation (Figure 4). 

Figure 4: improved model of LTD induction includes RKIP binding to Raf and MEK. Molecule numbers are indicated (from Hepburn et al. 2017).

An unresolved matter was an additional interaction of RKIP with Mitogen-activated protein kinase kinase (MEK), the kinase activated by Raf that activates ERK. It is known that RKIP also binds to and inhibits MEK, but while the RKIP-MEK interaction had been confirmed in Purkinje cells (Yamamoto et al. 2012) it was unclear whether this could block LTD induction. We therefore contacted Dr. K. Yamamoto (Korea Institute of Science and Technology) and she agreed to perform blocking experiments that confirmed that dissociation of the RKIP-MEK complex is required for LTD induction, as reported in our joint publication (Hepburn et al. 2017).
Several additional improvements were made in the model (Figure 4, Hepburn et al. 2017). The calcium profile caused by calcium uncaging was modeled more accurately, allowing for a better match with the data reported in Tanaka et al (2007). We also improved the phosphorylation of AMPAR by PKC, this resulted in a more robust model that fluctuates less than the one reported in Antunes and De Schutter 2012. 

Finally, in preparation of future spatial models, we explored the stochastic variability of the number of PKC molecules in the cytosol of a spine due to diffusion. For a mean of 55 PKC molecules (Figure 4), the predicted range is 40 to 70 molecules. Such changes are sufficient to shift the probabilistic calcium threshold for LTD induction by a factor of two (Figure 5). Therefore, diffusional mobility of key enzymes is an extra source of stochasticity in this system.

3.3 Information processing in the olivocerebellar system 

Multiplexed coding in primate Purkinje cells

The observation that PC simple spike trains in vivo consist of periods of regular firing interrupted by pauses has been suggested to represent a mixture of rate code (regular firing) with a temporal code (pauses). This would be an example of multiplexed coding, but the rodent data on which this hypothesis was based did not allow for behavioral correlation. Therefore, we explored whether we could confirm the multiplexed coding hypothesis in monkey data, which we obtained in a collaboration with Dr. Peter Thier (University of Tübingen, Germany).

We analyzed spike trains and local field potentials (LFPs) from 34 Purkinje cells obtained in three awake rhesus monkeys during spontaneous and visually guided saccades (Hong et al. 2016). We used an ISI asymmetry index (AI) instead of the CV2 because it distinguishes pause initiating from pause terminating spikes (Figure 6).  We confirmed that in all these PCs the spike trains could be described as a combination of prolonged periods of regular fast firing, occasionally interrupted by pauses (Figure 6), though the distinction was less sharp than in rodents. 

Figure 6: regular and pause spikes in the PC spike train. A: ISI asymmetry index (AI) measures local variability at each spike and is related to the local coefficient of variation (CV₂). B: distribution of AI for PC spikes (black) and rate-matching control spike train (green, mean ± SD). Significantly more PC spikes occurred around AI ≈ 0 (regular firing). C: three types of spiking pattern classified by their AI and associated ISI length; bottom: different types of spikes in an actual PC spike train (adapted from Hong et al. 2016).

We discovered that pause related spikes modulated differently with the LFP, this was shown by analyzing the Spike Triggered Average (STA) of the LFP. Pause-initiating and pause-terminating spikes showed strong phase-coupling to the LFP (Figure 7A), specifically to its β/γ band (Figure 7C), while regular spikes modulated only weakly with the LFP. This suggests that pause spikes are more strongly and more precisely coupled to network activity than regular spikes, but this fact would go unnoticed if the temporal structure of the spike train is ignored when computing the STA_LFP (Figure 7B), because pauses are relatively rare events. 

Figure 32: pause and regular spikes correlate differently with the LFP. A: STA_LFP of the different types of spikes. Gray region is the 99% confidence interval. B: STA_LFP of randomly selected spikes (black) and all spikes (green). Data are mean ± SEM. C: Phase distribution of the β/γ LFP at the pause, regular, and randomly selected spikes (thick: histogram, thin: kernel-estimated density). D: Pair phase consistency for each spike category (adapted from Hong et al. 2016).

If pause spikes in single PCs are related to specific patterns of network activity, what information do they jointly encode about eye motion? In many of PCs, the β/γ LFP was significantly coherent across all saccades, with a phase locking to saccade onsets (Figure 8A,B). This predicted that pause spikes, which are phase-locked to the β/γ LFP, should reliably code eye movement timing. The most significant pattern was that pause-initiating spikes encode saccade onsets with a significant and sharp increase in average firing, caused by an increase in spike coincidence (colored regions in Figure 8C). Spike reliability was significantly higher for pause spikes (Figure 8D), further emphasizing their role in encoding the onset of eye motion. Finally, the LFP phase and spike reliability were much more steeply related to each other for pause-initiating spikes than for regular spikes (Figure 8E). This demonstrates that the fidelity of temporal coding by pauses in individual PCs crucially reflects the reliability of activity and coding of the local network. 

Figure 33: pause spikes and β/γ LFP encode motion timing. A: eye speed and β/γ LFP aligned with saccade onset. Black lines represent an average over all saccades (gray). B: phases of the β/γ LFP at saccade onset. C Top: Spike trains for different types of spikes aligned with onsets of randomly selected saccades. Light colored regions represent periods with significantly reliable firing. Bottom: smoothed Z-score for spike occurrence. D: spike reliability of pause and regular spikes in all data. E: LFP phase reliability versus spike reliability (Hong et al. 2016).

So far we have focused on pause spikes, rare events in the spike train. On the other hand, the eye velocity-spike correlation observed in the data suggests that a firing rate code may also be present in PC spike trains. To answer this question, we constructed an inverse model that predicts firing rate from eye movements in each cell, based on the linear-nonlinear model framework. We found that the linear rate prediction followed the actual rate very closely. Importantly, correlation with the firing rate remained high when we computed the linear prediction with only regular spikes, which comprised about 20% of all spikes (Figure 9A). Conversely, the firing rate of pause spikes did not modulate as linearly or steeply as for regular spikes. This suggests that information encoded by the whole firing rate is very different from eye movement features related to pause spikes such as onsets (Figure 8C), but this information can be captured well by a subset of the full spike train, regular spikes (Figure 9B-C). 

Figure 9: PC firing rate linearly encodes motion kinematics as predicted by a linear model. A: predicted vs. measured firing rate for all (green), regular (black), and pause (red) spikes. The rates of pause and regular spikes are rescaled to match those of all spikes to compensate for subsampling. The dotted line represents equality. B: goodness of fit R₂ for the linear prediction to actual rate. C: R₂ for all cells (adapted from Hong et al. 2016).

The conclusion of our study is that PC spike trains can simultaneously convey information necessary to achieve precision in timing, as a temporal code linked to pauses in firing, and in continuous control of motion, as a rate code formed by regular spiking.

4. Publications

4.1 Journals

1. Chen, W. & De Schutter, E.  Parallel STEPS: Large Scale Stochastic Spatial Reaction-Diffusion Simulation with High Performance Computers. Frontiers in Neuroinformatics 11, doi:http://dx.doi.org/10.3389/fninf.2017.00013 (2017).

2. Chen, W. & De Schutter, E.  Time to Bring Single Neuron Modeling into 3D. Neuroinformatics 15, 1-3, doi:http://dx.doi.org/10.1007/s12021-016-9321-x (2017).

3. Gallimore, A. R., Aricescu, A. R., Yuzaki, M. & Calinescu, R.  A Computational Model for the AMPA Receptor Phosphorylation Master Switch Regulating Cerebellar Long-Term Depression. PLoS Computational Biology 12, e1004664, doi:10.1371/journal.pcbi.1004664 (2016).

4. Gallimore, A. R. & Strassman, R.  A Model for the Application of Target-Controlled Intravenous Infusion for a Prolonged Immersive DMT Psychedelic Experience. Frontiers in Pharmacology 7, 211, doi:10.3389/fphar.2016.00211 (2016).

5. Hepburn, I., Chen, W. & De Schutter, E.  Accurate reaction-diffusion operator splitting on tetrahedral meshes for parallel stochastic molecular simulations. Journal of Chemical Physics 145, 054118, doi:10.1063/1.4960034 (2016).

6. Hepburn, I., Jain, A., Gangal, H., Yamamoto, Y., Tanaka-Yamamoto, K. & De Schutter, E.  A Model of Induction of Cerebellar Long-Term Depression Including RKIP Inactivation of Raf and MEK. Frontiers in Molecular Neuroscience 10, 19, doi:http://dx.doi.org/10.3389/fnmol.2017.00019 (2017).

7. Hong, S., Negrello, M., Junker, M., Smilgin, A., Thier, P. & De Schutter, E.  Multiplexed coding by cerebellar Purkinje neurons. eLife 5, doi:10.7554/eLife.13810 (2016).

8. Negrello, M. & De Schutter, E.  Models of the Cortico-cerebellar System. Neuroscience in the 21st Century, 1-24 (2016).

4.2 Books and other one-time publications

Hepburn, I.  Realistic molecular simulation of neuronal signalling systems PhD thesis (2016).

4.3 Oral and Poster Presentations

Oral Presentations

1. De Schutter, E.  Modeling stochastic processes in neurons, KAIST, Korea (2016).

2. De Schutter, E.  Accurate Reaction-Diffusion Operator Splitting on Tetrahedral Meshes for Parallel Stochastic Molecular Simulations, Isaac Newton Instittute for Mathematical Sciences, Cambridge, UK (2016).

3. De Schutter, E.  Multiplexed coding by Purkinje cell simple spikes, KIST, Korea (2016).

4. Hong, S.  GABA-mediated phase repulsive coupling and seasonal time coding in the SC, KIST, Korea (2016).

5. Hong, S.  Multiplexed coding by cerebellar Purkinje neurons, Dept. of Brain & Cognitive Sciences, DGIST (Daegu, South Korea) (2016).

6. Hong, S.  Multiplexed coding in the cerebellar cortex, Annual Meeting of Korean Society for Chemical Senses (Ansan, South Korea) (2016).

7. Hong, S. & Myung, J.  GABA-mediated phase couplings and seasonal time coding in the suprachiasmatic nucleus, Kyushu University, Fukuoka Japan (2016).

Conference Invited lectures

1. De Schutter, E.  Modeling stochastic processes in neurons, in Volga Neuroscience Meeting 2016, Saint Petersburg-Nizhny Novgorod, Russia (2016).

2. De Schutter, E.  Multi-scale modeling of the cerebellum, in Collaborative Development of Data-Driven Models of Neural Systems, Janelia Research Campus, USA (2016).

3. Hepburn, I. & Gallimore, A. R.  Subcellular modeling tutorial, in CNS2016 Meeting, Jeju, Korea (2016).

4. Hong, S. & Myung, J.  GABA-mediated phase couplings and seasonal time coding in the suprachiasmatic nucleus, in JSMB2016 Fukuoka, Fukuoka, Japan (2016).

Poster Presentations

1. Chen, W. & De Schutter, E.  Toolkit Support for Complex Parallel Spatial Stochastic Reaction-Diffusion Simulation in STEPS, in CNS2016 Meeting, Jeju, Korea (2016).

2. Gallimore, A. R.  A unified molecular mechanism of bidirectional plasticity at cerebellar parallel fiberPurkinje cell synapses, Neuroscience2016 Meeting, San Diego, USA (2016).

3. Hong, S., Takashima, A. & De Schutter, E.  Dendritic morphology determines how dendrites are organized into functional subunits, in CNS2016 Meeting, Jeju, Korea (2016).

4. Hong, S., Takashima, A. & De Schutter, E.  Impact of dendritic morphology on functional subunits in dendrites,  Neuroscience2016 Meeting, San Diego, USA (2016).

5. Klingbeil, G. & De Schutter, E.  Stochastic Engine for Pathway Simulations on massivley parallel processors, in CNS2016 Meeting, Jeju, Korea (2016).

6. Torben-Nielsen, B., Bas, E., Chen, W., Cuntz, H., Kim, J., Kubota, Y., Moore, A. W., Shih, C.-T., Tavosanis, G., Peng, H., Ascoli, G. A. & De Schutter, E.  Future-proof digital representation of neuronal morphologies,  Neuroscience2016 Meeting, San Diego, USA (2016).

7. Verduzco Flores, S. & De Schutter, E.  Design of biologically-realistic simulations for motor control, in CNS2016 Meeting, Jeju, Korea (2016).

8. Zamore, C., Gallimore, A. R. & De Schutter, E.  A model of Ca2+/calmodulin-dependent protein kinase II activity in Long Term Depression at Purkinje cells, in CNS2016 Meeting, Jeju, Korea (2016).

9. Zang, Y. & De Schutter, E.  Modeling the generation and propagation of Purkinje cell dendritic spikes caused by parallel fiber synaptic input in CNS2016 Meeting, Jeju, Korea (2016).

 

5. Intellectual Property Rights and Other Specific Achievements

Nothing to report

6. Meetings and Events

6.1 OIST Mini-Sympo

Digital representation of neuronal morphologies and tissue

• Date: April 11-12, 2016
• Venue: OIST Campus C210, Center Building
• Organizer: Erik De Schutter (OIST)
• Speakers: 
  • Giorgio Ascoli (George Mason University, USA)
  • Erhan Bas (Janelia Research Campus, USA)
  • Weiliang Chen (OIST)
  • Hermann Cuntz (Ernst Strüngmann Institute (ESI), Frankfurt, Germany)
  • Jinny Kim (KIST, Korea)
  • Yoshiyuki Kubota (NIPS, Japan)
  • Adrian W. Moore (RIKEN BSI, Japan)
  • Hanchuan Peng (Allen Brain Institute, USA)
  • Chi-Tin Shih (National Tsing Hua University, Taiwan)
  • Gaia Tavosanis (German Center for Neurodegenerative Diseases, Germany)
  • Ben Torben-Nielsen (University of Hertfordshire, UK)

6.2 Summer Course

OIST Computational Neuroscience Course 2016

• Date: June 13-30, 2016
• Venue: OIST Seaside House
• Organizers: Erik De Schutter (OIST), Kenji Doya (OIST), Bernd Kuhn (OIST), and Jeff Wickens (OIST)
• Speakers:
  • Claudia Clopath (Imperial College London)
  • Erik De Schutter (OIST)
  • Kenji Doya (OIST)
  • Chris Eliasmith (University of Waterloo, Canada)
  • Tomoki Fukai (RIKEN BSI, Japan)
  • Michael Häusser (University College London, UK)
  • Etienne Koechlin (École Normale Supérieure, France)
  • Bernd Kuhn (OIST)
  • Stefan Mihalas (Allen Institute for Brain Science, USA)
  • Partha Mitra (Cold Spring Harbor, USA)
  • Shinji Nishimoto (CiNet, Japan)
  • Astrid Prinz (Emory University, USA)
  • John Rinzel (New York University, USA)
  • Greg Stephens (OIST)
  • Yoko Yazaki-Sugiyama (OIST)

6.3 Seminars

Title: Astrocytic GABA: biosynthesis, release, and function

• Date: January 18, 2017
• Venue: D015, Lab1, OIST
• Speaker: Prof. Changjoon Justin Lee, KIST (Korea Institute of Science and Technology)

 

Title: Computational modeling of astroglia-neuron interactions

• Date: March 09, 2017
• Venue: D014, Lab1, OIST
• Speaker: Prof. Marja-Leena Linne, Tampere University of Technology, Finland

7. Other

Nothing to report.