ORIGINAL RESEARCH article
Reconstruction and Simulation of a Scaffold Model of the Cerebellar Network
- 1Neurophysiology Unit, Neurocomputational Laboratory, Department of Brain and Behavioral Sciences, University of Pavia, Pavia, Italy
- 2IRCCS Mondino Foundation, Pavia, Italy
Reconstructing neuronal microcircuits through computational models is fundamental to simulate local neuronal dynamics. Here a scaffold model of the cerebellum has been developed in order to flexibly place neurons in space, connect them synaptically, and endow neurons and synapses with biologically-grounded mechanisms. The scaffold model can keep neuronal morphology separated from network connectivity, which can in turn be obtained from convergence/divergence ratios and axonal/dendritic field 3D geometries. We first tested the scaffold on the cerebellar microcircuit, which presents a challenging 3D organization, at the same time providing appropriate datasets to validate emerging network behaviors. The scaffold was designed to integrate the cerebellar cortex with deep cerebellar nuclei (DCN), including different neuronal types: Golgi cells, granule cells, Purkinje cells, stellate cells, basket cells, and DCN principal cells. Mossy fiber inputs were conveyed through the glomeruli. An anisotropic volume (0.077 mm3) of mouse cerebellum was reconstructed, in which point-neuron models were tuned toward the specific discharge properties of neurons and were connected by exponentially decaying excitatory and inhibitory synapses. Simulations using both pyNEST and pyNEURON showed the emergence of organized spatio-temporal patterns of neuronal activity similar to those revealed experimentally in response to background noise and burst stimulation of mossy fiber bundles. Different configurations of granular and molecular layer connectivity consistently modified neuronal activation patterns, revealing the importance of structural constraints for cerebellar network functioning. The scaffold provided thus an effective workflow accounting for the complex architecture of the cerebellar network. In principle, the scaffold can incorporate cellular mechanisms at multiple levels of detail and be tuned to test different structural and functional hypotheses. A future implementation using detailed 3D multi-compartment neuron models and dynamic synapses will be needed to investigate the impact of single neuron properties on network computation.
The causal relationship between components of the nervous system at different spatio-temporal scales, from subcellular mechanisms to behavior, still needs to be disclosed, and this represents one of the main challenges of modern neuroscience. The issue can be faced using bottom-up modeling, which allows propagating microscopic phenomena into large-scale networks (Markram, 2012; Markram et al., 2015; D'Angelo and Wheeler-Kingshott, 2017). This reverse engineering approach integrates available details about neuronal properties and synaptic connectivity into realistic computational models and allows to monitor the impact of microscopic variables on the integrated system. Realistic modeling can predict emerging collective behaviors producing testable hypotheses for experimental and theoretical investigations (Llinas, 2014) and might also play a critical role in understanding neurological disorders (Soltesz and Staley, 2018). In practice, realistic modeling of microcircuit dynamics and causal relationships among multi-scale mechanisms poses complex computational problems. First, the modeling strategy needs to be flexible accounting for a variety of neuronal features and network architectures, to be easy to update as new anatomical, or neurophysiological data become available, and to be easy to modify in order to test different structural and functional hypotheses. Secondly, the modeling tools need to be scalable to the dimension of the network and to the nature of the scientific question (Destexhe et al., 1996), to be suitable for available simulation platforms, e.g., pyNEST and pyNEURON (Brette et al., 2007; Eppler et al., 2008; Hines et al., 2009), and to efficiently exploit High-Performance Computing (HPC) resources.
Markram et al. recently carried out a digital reconstruction of the neocortical microcolumn by integrating experimental measurements of neuronal morphologies, layer heights, neuronal densities, ratios of excitatory to inhibitory neurons, morphological and electro-morphological composition, electrophysiological properties of neurons, and synapses (Markram et al., 2015). Neuron parameters were derived from databases specifically addressing cerebro-cortical neuron properties (e.g., Blue Brain Project and Allen Brain Atlas; Markram, 2006; Sunkin et al., 2013). Microscopic network wiring was then estimated computationally through a touch detection algorithm, that is based on a probability/proximity rule (i.e., the probability that morphologically defined dendrites and axons make a synaptic connection depends on their spatial proximity). This approach, in which the reconstruction of microcircuit connectivity depends on the 3D morphology of the axonal and dendritic processes of individual cells, may apply to brain structures for which datasets comparable to neocortex are available. However, such specific datasets are not available in general for all brain regions and it seems convenient in principle to keep separated neuronal morphology from network connectivity, which is reported as convergence/divergence ratios and axonal/dendritic field geometries in the literature in many cases.
The cerebellum hosts the second largest cortical structure of the brain and contains about half of all brain neurons. Modeling the cerebellum brings about specific issues reflecting the peculiar properties of this circuit, which shows a quasi-crystalline geometrical organization well-defined by convergence/divergence ratios of neuronal connections and by the anisotropic 3D orientation of dendritic and axonal processes (Figure 1) (D'Angelo et al., 2016). Moreover, the morphological reconstruction of axonal and dendritic processes of cerebellar neurons is not as developed as for other brain microcircuits, like cerebral cortex, and hippocampus (e.g., see the NeuroMorpho.org repository —https://www.re3data.org/repository/r3d100010107; Akram et al., 2018). Therefore modeling the cerebellum relies on a knowledge base that differs from that available for the cerebral cortex and thus requires a more general approach than in the Markram's modeling workflow (Markram et al., 2015).
Figure 1. Reconstruction of a scaffold model of the cerebellar network. Schematic representation of the cerebellar network (from D'Angelo et al., 2016). Glomeruli (Glom); mossy fiber (mf); Granule cells (GrC); ascending axon (aa); parallel fiber (pf); Golgi cells (GoC); stellate cell (SC); basket cell (BC); Purkinje cell (PC); Deep Cerebellar Nuclei cell (DCNC). Gloms transmit mf inputs to GrCs, which emit aa and pf, which in turn activate GoCs, PCs, SCs, and BCs. GoCs inhibit GrCs, SCs and BCs inhibit PCs. DCN cells are inhibited by PCs and activated by mf. Note the precise organization of PC dendrites, SC/BC dendrites and GoC dendritic arborization mainly on the parasagittal plane. The same abbreviations are used also in the following figures.
Some recent models were purposefully designed to reproduce a limited section of the cerebellar cortex, the granular layer (Maex and De Schutter, 1998; Solinas et al., 2010; Sudhakar et al., 2017), in great detail and incorporated Hodgkin-Huxley-style mechanisms and neurotransmission dynamics (D'Angelo et al., 2001; Solinas et al., 2007a,b; Nieus et al., 2014; Masoli et al., 2015, 2017; Masoli and D'Angelo, 2017). Other models were designed to simulate, in a simplified form, large-scale computationally efficient networks of the olivo-cerebellar system (Medina and Mauk, 2000; Yamazaki and Nagao, 2012). In this work, a new cerebellum scaffold model has been developed and tested, allowing to incorporate axonal/dendritic field geometries specifically oriented in a 3D space and to reconnect neurons according to convergence/divergence ratios typically well-defined for the cerebellum (D'Angelo et al., 2016). The cerebellum scaffold model maintains scalability and can be flexibly handled to incorporate neuronal properties on multiple scales of complexity and to change its connectivity rules. For the sake of simplicity, here we used first simplified neuron and synaptic models to evaluate the impact of construction rules. The cerebellum scaffold model was validated by testing its ability to reproduce the structural properties anticipated experimentally and the emergence of complex spatiotemporal patterns in network activity. The model was run on the pyNEST and pyNEURON simulators (Brette et al., 2007; Eppler et al., 2008; Hines et al., 2009) and a test workflow was integrated into a large-scale neuroinformatics infrastructure, the Brain Simulation Platform (https://collab.humanbrainproject.eu/).
Materials and Methods
This paper reports the design and implementation of a scaffold model of the cerebellum microcircuit. The model architecture is scalable and is designed to host different types of neuronal models and to determine their synaptic connectivity from convergence/divergence ratios and axonal/dendritic field geometries reported in literature. The workflow encompasses two main modules in cascade: cell placement into a user-defined volume; connectivity among neurons. The scaffold can then be used for functional simulations of network dynamics (Figure 1).The scaffold is designed to be embedded into different simulators, e.g., pyNEST and pyNEURON. This workflow, by allowing a flexible reconstruction of the cerebellar network, will eventually allow to evaluate physiological, and pathological hypotheses about circuit functioning.
Cell Placement Module
The Cell Placement Module placed the cells in a virtual network volume divided in layers based on morphological definitions. The process took into consideration the different cerebellar neuron types: the Golgi Cell (GoC), Granule Cell (GrC), Purkinje Cell (PC), Stellate Cell (SC), Basket Cell (BC), Deep Cerebellar Nuclei glutamatergic GAD-positive Cell (DCNC), and glomerulus (Glom). Glom is actually a mossy fiber terminal and is represented as a neuronal element at the input stage, while DCNCs are placed at the output stage of the circuit. For each neuron type, the density value in a specific layer was derived from literature, and geometric features (including soma radius and 3D-oriented dendritic and axonal fields) were defined according to experimental data. Through ad hoc algorithms (Bounded Self-Avoiding Random Walk Algorithm and Purkinje Cells placement algorithm, see below; and details in Supplementary Material Placement workflow), the cells were positioned in the 3D volume of each layer, according to their density, soma radius, and anisotropic extension, ensuring that their somata did not overlap. The module was implemented in Python, and its output was saved in an. hdf5 file containing the unique identification number (ID) of each cell, its corresponding type (an integer value between 1 and 7, as in Table 1), and the three spatial coordinates of the soma center (x, y, z). To evaluate the effectiveness of the cell positioning algorithms, we derived a continuous distribution of pair-wise distances using kernel density estimation (KDE), in which the Gaussian kernel had fixed bandwidth for each cell population. KDE yielded a single maximum when pair-wise distances were distributed homogeneously (GrC, GoC, SC, BC, DCNC) and multiple local maxima when distances were placed according to different geometric rules (PC). A reconstructed network volume and pair-wise soma distances yielded by this module are illustrated in Figure 2.
Figure 2. Cell placement and network architecture. (A) The cells are placed in the network 3D space using a Bounded Self-Avoiding Random Walk Algorithm. The figure shows the volume of 400 × 400 × 900 μm3 containing 96,737 neurons and 4,220,752 synapses used for simulations. (B) Projection of GrC axons to the molecular layer hosting the PCs (green dots in the PC layer are the somata, the thin green parallelepipeds above are the corresponding dendritic trees occupying the molecular layer). The figure shows two clusters of GrCs and the corresponding aa and pf, illustrating that the cerebellar network connectivity respects the 3D architecture shown in Figure 1. (C) Distributions of 3D pair-wise inter-soma distances within each neuronal population: GrCs, SCs, GoCs, BCs, and PCs. Note that the distributions are nearly normal, except for the PCs.
GrCs, GoCs, SCs, BCs, and DCNCs were placed in thin sublayers (with height = 1.5x soma diameter) using a bounded self-avoiding random walk algorithm. In each sublayer, the cells were initially distributed in 2D and then sublayers were piled one on top of the others. The first cell was placed randomly and each subsequent one was positioned nearby along a random angular direction. The overlap among somata was avoided since, along the selected direction, the minimum distance to place the next cell was equal to the soma diameter. A term was added to the minimum distance to scatter the somata. In details, a potential volume for each cell was computed from density values, then deriving the difference between this compound sphere radius and the soma radius (ε); a value was randomly sampled from as a normal distribution around ε (minimum 0.75· ε, maximum 1.25· ε). This guarantees natural randomness but at the same time a good exploitation of the whole available volume, avoiding undesired clusters or not uniform occupancy. If the surrounding space was completely occupied, i.e., there was insufficient space to place a further cell, a new starting point was selected resetting the random walk process for the remaining neurons in that sub-layer. Once completed, the 2D sub-layer was piled on top of the underlying one. Then, a vertical coordinate was assigned to each cell, from a random uniform distribution within the sublayer height (Korbo et al., 1993). This approach maintained randomness achieving a realistic quasi-Gaussian distribution of pair-wise inter-neuron distances (see Figure 2C) and proved computationally efficient.
The PCs were distributed in a single sub-layer forming an almost planar grid between the granular and molecular layers. The PC inter-soma distances over this plane were constrained by the dendritic trees, which are flat and expand vertically on the parasagittal plane (about 150 μm radius × 30 μm width) without overlapping (Masoli et al., 2015). Since PC somata do not arrange in parallel arrays but are somehow scattered, a noisy offset was introduced creating an average angular shift of about 5° between adjacent PCs. As for the other neuron types, a small random noise was also imposed on the vertical coordinate (Korbo et al., 1993).
The data required for cell positioning in the cerebellar cortex were obtained from literature (Eccles et al., 1967; Magyar et al., 1972; Mezey et al., 1977; Hamori and Somogyi, 1983; Jakab and Hamori, 1988; Korbo et al., 1993; Sultan, 2001; Santamaria et al., 2007; Barmack and Yakhnitsa, 2008; Solinas et al., 2010) and are summarized in Table 1. GoCs, GrCs, and Gloms were placed into the granular layer; BCs and SCs in the lower and upper half of the molecular layer, respectively. A certain number of DCNCs was randomly distributed in DCN volume according to the PC/DCNC ratio, since more specific parameters are still missing (Gauck and Jaeger, 2000; Aizenman et al., 2003; Person and Raman, 2012). Special care was given to the GrC ascending axon (aa) that, starting directly from the soma, makes its way up vertically toward the molecular layer. The height of each ascending axon was chosen from a Gaussian distribution in the range of 181 ± 66 μm (Huang et al., 2006). This value represents also the vertical coordinate of the corresponding parallel fiber (pf), running transversally and parallel to the cerebellar surface.
The connectivity module created structural connections between pairs of neurons belonging to specific types. Each neuron type formed input and output connections with other neurons of the same or different types. Therefore, once the placement was completed, it was possible to reconstruct the connectome applying intersection-connectivity rules based on proximity of neuronal processes and on statistical ratios of convergence and divergence. When available, morphological and statistical literature data were used, otherwise plausible physiological constraints were applied. In our scaffold, 16 connection types were generated (the most important are shown in Figure 3), from the volume covered by pre-synaptic axonal processes to that covered by post-synaptic dendritic trees of specific neuron types:
1. From glomeruli to granule cells (Glom-GrC);
2. From glomeruli to basolateral dendrites of Golgi cells (Glom-GoC);
3. From Golgi cells to Gloms (GoC-Glom): this is fused together with Glom-GrC to generate directly GoC-GrC connections;
4. Among Golgi cells (GoC-GoC);
5. From ascending axons of granule cells to Golgi cells (aa-GoC);
6. From parallel fibers of granule cells to apical dendrites of Golgi cells (pf-GoC);
7. Among stellate cells (SC-SC);
8. Among basket cells (BC-BC);
9. From parallel fibers of granule cells to stellate cells (pf-SC);
10. From parallel fibers of granule cells to basket cells (pf-BC);
11. From stellate cells to Purkinje cells (SC-PC);
12. From basket cells to Purkinje cells (BC-PC);
13. From ascending axons of granule cells to Purkinje cells (aa-PC);
14. From parallel fibers of granule cells to Purkinje cells (pf-PC);
15. From Purkinje cells to DCN cells (PC-DCNC);
16. From glomeruli to DCN cells (Glom-DCNC).
Figure 3. Cell connectivity: examples for specific connections. Examples of divergence and convergence at different connections in the cerebellar network space. The plots have base area (400 × 400 μm2) and thickness specific for each layer. The plots show a randomly selected pre-synaptic cell together with its connected post-synaptic neurons (divergence) or viceversa (convergence). (A) Connections of GrCs and GoCs. (B) Connections of PCs, SCs, and BCs. (C) Connections of DCNCs.
Given a connection type, for each pre-synaptic neuron, the potential post-synaptic cells were identified as those that met geometric neuron-specific constraints. Then, given the convergence and divergence ratios, post-synaptic neurons were selected among the potential ones, through a pruning process using distance-based probability functions specific for each volume direction (details and examples in Figure S1). The module was implemented in Python, and its output saved in an. hdf5 file containing a matrix for each connection type, in which each row was defined by three values: the unique ID of the pre-synaptic neuron, the unique ID of the post-synaptic neuron and the inter-soma 3D distance between that pair (see Figure 4A).The plots in Figures 4B,C compare, for each connection type, the divergence and convergence ratios reported by literature to the values obtained after scaffold reconstruction in a sample volume. The cell placement and connections rules yielded indeed a very good approximation of the anatomical and physiological parameters reported in literature.
Figure 4. Cell connectivity: pair-wise distance prediction and convergence/divergence validation. (A) Pair-wise distance prediction deriving from the placement and subsequent cell-to-cell connectivity. The data that find correspondence in literature are indicated as asterisks. For each connection type, the pair-wise distances between connected cells (inter-soma distance) are reported. Data from: (1) (D'Angelo et al., 2013), (2) (Barmack and Yakhnitsa, 2008), (3) (Rieubland et al., 2014). (B,C) The plots compare divergence and convergence for the different connections of the scaffold with those anticipated experimentally. The regression lines show a very close correspondence of the model to experimental results. Linear regression lines are fitted to the data (divergence: r2 = 0.98, slope = 0.88; convergence: r2 = 0.99, slope = 0.99). Data from: (1) (Nieus et al., 2006), (2) (Dieudonne, 1998), (3) (D'Angelo et al., 2013), (4) (Solinas et al., 2010), (5) (Kanichay and Silver, 2008), (6) (Cesana et al., 2013), (7) (Hull and Regehr, 2012), (8) (Lennon et al., 2014), (9) (Huang et al., 2006), (10) (Jorntell et al., 2010), (11) (Sultan and Heck, 2003), (12) (Person and Raman, 2012), (13) (Boele et al., 2013).
In order to test the functionality of the scaffold, single neuron models were placed in the corresponding positions of the connectome matrix. In this first version of the cerebellar microcircuit, spiking point-neuron models based on Integrate&Fire (I&F) dynamics with conductance-based exponential synapses (i.e., synaptic inputs cause an exponential-shaped change in synaptic conductances) were used. The output files of these simulations contained all the spike events (neuron IDs and relative spike times). Glomeruli were represented as “parrot neurons” just able to pass the imposed stimulation patterns unaltered. Each other neuron type was characterized by specific values, directly related to neurophysiological quantities (Cm, τm, EL, Δtref, Ie, Vr, Vth), corresponding to biological values taken from literature available from animal experiments or databases (https://neuroelectro.org/) (Tripathy et al., 2014; Table 2). In order to account for the neuron-specific dynamics of GABA and AMPA receptor currents, also the decay times of the excitatory and inhibitory synaptic exponential functions (τexc, τinh) were set differently for each neuron type (Table 2). Each synaptic connection type was characterized by specific values of weight and delay (Table 3). These estimates approximated literature data values so that, for example, the synaptic delay was shorter when fibers impinging on PCs came from aa than pf synapse.
The input stimulus was set by defining the volume where Gloms were activated, the onset time, the duration and the frequency of spikes. A background activity was generated by a Poisson process of stochastic neuronal firing at 1 Hz on all the glomeruli. Superimposed on it, a burst at 150 Hz lasting 50 ms (Rancz et al., 2007) was activated 300 ms after the beginning of simulation. Indeed, mossy fibers have a low basal activity, but in response to sensorimotor stimuli, can fire at rates beyond 100 Hz. The stimulated volume had a radius of 140 μm; the simulation lasted 1 sec, including 300 ms pre-stimulus, 50 ms stimulus, and 650 ms post-stimulus.
In a specific set of simulations (see Figure 8), we tested the partial contribution of SCs and BCs to the spatiotemporal diffusion of activity among PCs. BCs axonal plexus is preferentially oriented along the parasagittal axis (see Eccles et al., 1967). In these simulations, we oriented the SC and BC axonal plexus orthogonally one to each other and concentrated the stimulation burst in a sphere of 30 μm radius.
Network Data Analysis
For each neuron population, mean frequency rates were extracted in three time windows: baseline pre-stimulus, during stimulus, and steady-state after-stimulus. We then generated peri-stimulus time histograms (PSTH) for each neuronal population with time bins of 3 ms. For each neuron population, we also separated excited from inhibited sub-groups, responding with an increased, or decreased firing rate during the stimulus. To do so, we compared the number of spikes during stimulus vs. pre-stimulus normalized by the time-window durations. If the pre-stimulus firing frequency (i.e., baseline) was at least doubled during stimulation, then the cell was classified as excited. Conversely, to classify the inhibited cells. For GrCs, we added a second constraint: to be labeled as excited, a GrC should fire more than 1 spike during stimulation, allowing to exclude spikes determined by the background noise. For DCNC, all cells stopped firing during the stimulation time-window. For each PC, a further ad-hoc analysis allowed to identify burst–pause responses. The cells showing a significant stimulus-induced pause (Cao et al., 2012; Herzfeld et al., 2015) were recognized as those in which the first Inter-Spike-Interval (ISI) after the end of the stimulus was >2 standard deviation (sd) of the pre-stimulus ISIs. This comparison was computed within-cell, i.e., for each PC individually.
The Excitatory-Inhibitory balance (EI) and Center-Surround (CS) were calculated from firing rates (FR) according to Mapelli and D'Angelo (2007) and Solinas et al. (2010) by considering that inhibition occurs only after a delay following the beginning of stimulation. GrC firing rate was then measured 0–20 ms (T1) and 20–40 ms (T2) after the beginning of stimulation in response to 50 ms at 150 Hz bursts, both in control (con)and with GoC-GrC inhibition switched off (in_off ). The CS and EI were calculated as follows:
The CS values were normalized between 1 and −1. The extension of the center and surround was calculated by including zones with CS > 0.5 in the center, and zones with CS < −0.5 in the surround. The center and surround relative areas could then be calculated by counting the respective number of pixels and normalizing by the total number of pixels (see Figure 7C).
In order to determine the presence and properties of coherent oscillations in granular layer activity, the activity in a subset of GoCs with overlapping incoming parallel fibers and the related GrCs was analyzed (Maex et al., 2000) during a 5 s at 5 Hz noisy background mossy fiber activity. The autocorrelations of GoCs and GrCs spike trains and the cross correlation between GrCs and GoCs spike trains were calculated using the equation
Where C is the index of coherence, A is the array of autocorrelation values, and len(A) is the size of the spike train data array. The same calculus was executed also for the crosscorrelation.
Simulations in pyNEST, pyNEURON, and Implementation on the Brain Simulation Platform
The microcircuit was implemented and simulated both in pyNEST version 2.14 (Eppler et al., 2008) and in pyNEURON (Hines et al., 2009). These neural simulators are commonly used for applications starting from realistic neuron models and up to more abstract representations. These tests were run using external HPC and local resources, maximizing available parallel computing. The time resolution for both simulators was 0.1 ms. As internal validation tests, some exemplificative ad-hoc structural and functional alternatives (see Figure 5) were made in the network and then the same simulations were run (in pyNEST). The firing rates of each cell population and their sub-groups affected by stimulus are reported in Table 4 and Figure 5 to illustrate the spiking network behaviors.
Figure 5. Neuronal discharge. Raster plot and PSTH of the different neuron populations of the cerebellar network model in response to a mossy fiber burst (50 ms at 150 Hz on 2,932 gloms) superimposed on a 1 Hz random background. The two simulations used the same cerebellar scaffold and neurons, which were translated from pyNEST into pyNEURON. The basal activity of the different cell populations is visible before and after the stimulus. The Glom patterns at the input are imposed, so they are identical for both simulations. The mean population firing rates for GrCs are similar between the two simulations, probably due to the very high number of GrCs. Minor differences are detectable for the other neuron types.
The entire scaffold can be built and run as a Jupyter Notebook in the Brain Simulation Platform (BSP), one of the platforms of Human Brain Project (Markram, 2012). The BSP is an internet-accessible collaborative platform that comprises a suite of software tools and workflows to reconstruct and simulate multi-level models of the brain at different levels of description, from abstract to highly detailed. Here, cells, network, and volume configuration parameters can be easily read and modified, since they are stored in a single Python script. Such flexible parametric approach allows to continuously include and tune relevant neurophysiological information and to operate at different simplification levels. A test version of the scaffold model is running on the Brain Simulation Platform at https://www.humanbrainproject.eu/en/brain-simulation/brain-simulation-platform/.
The cerebellar network is unique for its precise geometrical organization (Figure 1), which was reconstructed generating a scaffold model capable of handling neuronal placement, connectivity, and simulations. The neurons were represented as single-point leaky integrate-and-fire (LIF) models (Maas, 1997), tuned to match the input resistance and capacitance, basal discharge, and input-output relationship of the specific cerebellar neuron types. The choice of LIF neuron models was motivated by the need to focus first onto the two main construction operations of the scaffold, cell placement and connectivity, and on the role of these latter in determining network properties.
The scaffold is demonstrated here through the exemplar reconstruction and testing of a cerebellar volume of 0.077 mm3. The cerebellar cortex volume had 400 × 400 μm2 base and 330 μm height subdivided into different layers: molecular layer (150 μm), Purkinje cell layer (30 μm), granular layer (150 μm). The DCN layer had 200 × 200μm2 base (1/4 of cortex) and 600 μm height. As a whole, the model contained 96,734 cells. It should be noted that these parameters were all user-defined and may be modified depending on the needs, as the model is scalable.
The Bounded Self-Avoiding Random Walk algorithm (see section Materials and Methods) successfully placed the neurons into all cerebellar regions with the only exception of PCs, which were positioned using a specific algorithm designed to respect their regular spatial organization (Figure 2A). Figure 2B shows a row of almost equally distanced PCs connected to incoming parallel fibers, faithfully reproducing the typical PCs geometrical organization. These examples show that the placement algorithms can be flexibly configured to account for complex and variable rules of cellular positioning.
As an internal validation, the distribution of pair-wise distances for each cell type was calculated (Figure 2C). For all cell types (except PCs), pair-wise distances were distributed almost normally and the minimum inter-soma distance equated twice the soma radius. As expected, KDE for GrC, GoC, SC, and BC pair-wise distances returned a single maximum (at 180.1, 191.0, 184.6, and 188.5 μm, respectively), while for PCs three local maxima occurred (at 48.6, 142.1, and 267.0 μm) (for DCNC, KDE analysis was meaningless, given the low number of cells).
The connection rules adopted in this work were designed to account for the rich and specific information available from literature (Eccles et al., 1967; Palay and Chan-Palay, 1974; Korbo et al., 1993), which accounts for convergence/divergence ratios, number of synapses, and spatial distribution of axons and dendrites (Figure 3). The connecting algorithm imposed these geometrical constraints allowing to wire the different neuronal types for a whole of 16 connection types. Five connection types did not require other than these geometric constraints, while pruning was needed in the other 11 cases (either for convergence or divergence or both). The resulting connectome was then compared to the experimental one for validation. Figure 4 shows that the connection ratios of the scaffold were indeed correctly scaling to the physiological ones. Some specific cases are considered below.
Concerning Glom-GrC connectivity, experimental data demonstrated that granule cell dendrites have a maximum length of 40 μm, with a mean value of ~13 μm (Solinas et al., 2010). By imposing a convergence value of 4 (each GrC received one Glom on each of its 4–5 dendrites), a mean dendrite length of about ~12 μm was found, therefore matching experimental and theoretical determinations (Hamori and Somogyi, 1983; Billings et al., 2014).
Concerning connectivity between the aa and PC dendrites (aa-PC), connections were possible only when the aa-segment was very close to the PC dendritic plane. By analyzing the placement of GrCs in the x-z plane and the vertical extension of the aa, it is estimated that only ~20% of GrCs developed an aa that is sufficiently close to a PC dendrite tree to form a synaptic contact (Bower and Woolston, 1983; Gundappa-Sulur et al., 1999). This estimate was indeed closely matched by the scaffold reconstruction.
Concerning connectivity of parallel fibers with receiving neurons (pf-GoC, pf-SC, pf-BC, pf-PC synapses), the literature is incomplete and shows variable estimates. This most likely reflects difficulties in estimating exact numbers, since the pf can be several millimeters long and they are often cut on the parasagittal plane in histological preparations. In the scaffold reconstruction, the maximum pf length (along z-direction) was bounded to 400 μm (Barbour, 1993; Huang et al., 2006) and pf from GrCs beyond this length were not taken into account.
The statistical distribution of distances between connected cells (Table 4) shows a good matching with anatomical values. This validation of the connectome supports the appropriateness of cell placement and connecting rules. Biological randomness in the 3D placement with uniform occupancy of appropriate layers ensures that the resulting connectivity (based on geometrical proximity) has plausible biological values and variability for statistical convergence/divergence ratios, and for distances among connected neurons.
Neuronal Activations in the Cerebellar Network Following Mossy Fiber Stimulation
The aim of these simulations was to assess the emergence of typical spatio-temporal patterns of cerebellar network activity as a consequence of mossy fiber inputs. Simulations were carried out both in pyNEST and pyNEURON. For simplicity, the following data and figures are taken from pyNEST simulations, except for a comparison of the two in Figure 5 and Table 4. As expected, the two simulators yielded similar firing rates in each cell population. The Glom patterns at the input were imposed, so they were identical for both simulation platforms, while very small differences were detectable for the other neuron types.
Evoked activity simulating the effect of natural sensory stimulation (Chadderton et al., 2004; Roggeri et al., 2008; Ramakrishnan et al., 2016)was elicited over a noisy background (see above) by a 150 Hz−50 ms mossy fiber burst. The mossy fiber activity spread over about 0.012 μm3 of the granular layer involving 2,932 glomeruli out of the 7,070 placed in the reconstructed volume. Glomeruli had mean firing rate of ~1 Hz before the burst, 140 Hz during the burst, and ~1 Hz after the burst. The burst induced transient activity changes, specific for each neuronal population, which reverted back to baseline after the end of the stimulus (Figure 5). The sequence of neuronal activations depended on synaptic delays that were set according to physiological data (Eccles et al., 1967; Figure 6). The response of the individual neuronal populations was as follows (Figures 5, 6, and Table 4):
• The GrCs discharged at 1.8 Hz at rest and at 114 Hz during burst stimulation, consistent with in-vivo data showing that GrCs had sparse activity characterized by low background firing rates (partly due to the presence of tonic GABAergic inhibition) and high-frequency bursts in response to evoked sensory stimulation (Chadderton et al., 2004).
• GoCs discharged above 22 Hz at rest and above 150 Hz during burst stimulation, consistent with in-vivo data (Heine et al., 2010). The basal GoCs firing rate was raised by the noisy background over the autorhythmic frequency and showed a high variability among cells.
• Molecular layer interneurons, SCs and BCs (N = 603 for each cell type), discharged at ~30 Hz at rest and above 120 Hz during burst stimulation, consistent with the observation of high-frequency activity during sensory stimulation (Chu et al., 2012).
• PCs discharged at ~58–60 Hz at rest and at ~255 Hz during burst stimulation consistent with in-vivo data (Heine et al., 2010). Interestingly, PCs showed either bursts, or pauses, or burst-pause responses as observed in vivo (Herzfeld et al., 2015): out of 69 PCs, the burst was observed in 48 PCs and the pause in 41 PCs. Of the PCs that showed a pause, in 17 PCs it occurred after a burst, while in the other 24 PCs it happened alone.
• DCNCs discharged at ~16 Hz at rest and were completely silenced during burst stimulation. This behavior was expected from the convergent inhibition coming from PCs, supporting the hypothesis that cortico-nuclear synapses act as simplified inverters (Person and Raman, 2012).
Figure 6. Cerebellar network response to a mossy fiber burst. (A) Spikegrams of all cerebellar neurons in the model. A burst in gloms causes a burst-to-burst propagation in GrCs and PCs. GoCs, SCs, and BCs also generate bursts that, by being inhibitory, contribute to terminate the GrC, and PC bursts and to generate the burst-pause PC response. The DCN cells show a pause during stimulation. (B) Raster plot of one cerebellar neuron for each population in the model. Note the spread of the mf bursts inside the cerebellar cortical networks and the corresponding pause in the DCN. (C) Spike-time response plot showing the temporal sequence of neuronal activation and inhibition. The arrows represent the connectivity (solid lines show excitatory connections, dashed lines inhibitory connections). The stars represent the post-synaptic neuron response: white stars are excited neurons, black stars are inhibited neurons.
Center-Surround Organization of Granular Layer Responses
A relevant aspect of network activation that emerged in electrophysiological and imaging experiments is the center-surround organization (Mapelli and D'Angelo, 2007; Diwakar et al., 2011; Gandolfi et al., 2014). In the scaffold, the neuronal response of the granular layer to mossy fiber stimulation showed a typical center-surround organization (Figure 6). This reflected the excitatory/inhibitory ratio (see section Materials and Methods) with the center more excited than the surround due to lateral inhibition provided by GoCs. The center-surround had a diameter of about 50 μm and GrCs inside the core discharged up to 3–4 spikes organized in a short burst, reflecting previous experimental estimates (Gandolfi et al., 2014). Therefore, the scaffold correctly predicts the consequences of activity in bundles of mossy fibers.
Recently the connectivity of GoCs and GrCs has been extended by the demonstration of new synapses, in particular those between the GrC ascending axon and GoCs (aa-GoC, excitatory) (Cesana et al., 2013) and between GoCs (GoC-GoC, inhibitory) (Hull and Regehr, 2012). The selective switch-off of aa-GoC connections enhanced the center and reduced the surround, the switch-off of GoC-GoC connections reduced the center and increased the surround, while smaller effects followed the switch-off of pf-GoC or mf-GoC synapses (Figure 9C).
The Impact of Molecular Layer Interneurons on PC Activation
The molecular layer is critical to regulate PC activity in a way that is still debated (e.g., see Rokni et al., 2007; Santamaria et al., 2007). The first assumption is a differential orientation of SC cell axons (mostly transversal or “on-beam”) vs. BC axons (mostly sagittal or “off-beam”) (Eccles et al., 1967). Moreover, both aa and pf are used to activate PCs, as reported in literature (Jaeger and Bower, 1994; Canepari et al., 2001; Figure 8). Consistently, in the scaffold model, PC responses were circumscribed into a central spot overlaying the center/surround generated in the granular layer with little diffusion along either transversal or sagittal axis. On both axes, in turn, some PCs were clearly inhibited by the molecular layer interneuron inhibitory network.
Figure 7. Center-surround organization of activity in the granular layer. (A) In response to a mossy fiber burst (40 gloms at 150 Hz for 50 ms), the granular layer responds with a core (red area) of activity surrounded by inhibition (blue area). (B) PSTH of GrCs in the center-surround. The activity in the core is characterized by robust spike bursts, while just sporadic spikes are generated in the surround. No activity changes are observed outside the center-surround structure. (C) The histogram shows the changes in center-surround extension that occur following selective switch-off of synapses impinging on GoCs. Note the prominent role of aa-GoC synapses and GoC-GoC synapses (bars are values normalized to control).
Figure 8. Maps of PC activation and sensitivity to molecular layer connectivity. (A) The maps show the activity change of PCs in response to a mossy fiber burst (40 glom at 150 Hz for 50 ms). The pattern of activity is determined by various connection properties that are tested in turn. (all active) PC inhibition is achieved through a differential orientation of SC axons (mostly transversal or “on-beam”) vs. BC axons (mostly sagittal or “off-beam”) and that PC excitation depends on both aa and pf synapses with specific origin from GrCs. Alternative patterns are generated by (SC off) the specific switch-off of SC, (BC off) the specific switch-off of BC, or (SC&BC off) the complete switch-off of both SC and BC, (aa off) the specific switch-off of aa synapses, (pf off) the specific switch-off of pf synapses. It should be noted that these changes in network connectivity modify the PC discharge patterns both on-beam and off-beam and extend to a distance that reflects the propagation of activity through the molecular layer interneuron network. The circles indicate the location of the underlying active spots of activity in the granular layer. The bottom plot represents the activity of GoCs (blue) and GrCs (red) before, during and after the stimulus burst. This activity occurs in a spot (enlarged in the inset) corresponding to the center-surround shown in Figure 7. (B) The schematic diagrams show the orientation of fibers and connections in the network. (C) The PC activity was averaged into 3 × 3 matrices in order to better appreciate where activity changes take place. Note the emergence of the central spot in several cases.
Then, the effect of disconnecting different network elements was tested. Following the switch-off of both SC and BC inhibition, the responsiveness of PCs increased, as expected from SC and BC inhibitory action on PCs. As expected from anatomy, when only BCs were present (i.e., selective switch-off of SCs), excitation extended more effectively along the transverse axis, while when only SCs were present (i.e., selective switch-off of BCs) excitation extended more effectively along the sagittal axis. However, in both cases there was a diffused (though slight) increase of excitation, due to the reduced background inhibition exerted by intrinsic SC and BC discharge. It should also be noted that the activation of PCs in the central spot remained poorly altered, suggesting that these PCs were already nearly maximally activated in control. The selective switch-off of aa synapses caused a diffuse reduction of PC activation, while the selective switch-off of pf synapses had a much smaller effect. Therefore, changes in molecular layer connectivity consistently modified the PC discharge patterns both on-beam and off-beam and extended to a distance that reflects the propagation of activity through the pfs and the molecular layer interneuron network.
Synchronous Oscillations Caused by Noisy Background Activity in Mossy Fibers
Recordings from the granular layer in vivo have revealed low-frequency local field potential oscillations that occur synchronously over distances of several hinders of micrometers (Pellerin and Lamarre, 1997; Hartmann and Bower, 1998). Similar properties were observed also in previous granular layer models (Maex and De Schutter, 1998; Solinas et al., 2010; Sudhakar et al., 2017). In the scaffold model, spontaneous circuit activity clearly emerged due to background firing in the mossy fibers, provided that the frequency of the background mossy fiber discharge was increased from 1 to 5 Hz and pfs-GoCs connection weight was increased from 0.4 to 30.4, supporting the concept that oscillations require a specific synaptic balance to emerge (Maex and De Schutter, 1998; Solinas et al., 2010; Sudhakar et al., 2017; Figure 9). In response to the input, GrCs sparsely discharged at low frequencies (GrCs do not show intrinsic spontaneous activity), while the intrinsic activity of all the other neurons was modulated (GoC, PC, MLI, and DCNC are autorhythmic) (see Ie values in Table 3). Interestingly, the neurons of the granular layer (GrCs and GoCs) showed a pattern of low-frequency oscillations (mean frequency of 1.8 Hz) that was evident across the whole network. The oscillation frequency is the same in autocorrelograms of both GrCs and GoCs, and in the cross-correlogram between Golgi and granule cells. This ensemble behavior is probably due to the inhibitory feedback from GoCs to GrCs in the following way: (1) GrCs activity sums up in several GoCs, (2) GoCs, which are synchronized through parallel fibers and reciprocal inhibitory synapses, discharge almost synchronously, (3) a large population of GrCs is phasically inhibited, (4) inhibition terminates and GrCs recover responsiveness to the background mossy fiber input, then restarting the cycle.
Figure 9. Coherent low-frequency oscillations in granular layer neurons. Activity of GrCs (red) and GoCs (blue) during sustained 5 Hz random mf input. (A) Raster plots from exemplar GrCs and GoCs. Note that synchronous patterns are visible in the neuronal response (arrows). In this regimen, GoC activity is more intense than GrC activity due to the autorhytmic discharge of GoC neurons. The neurons are not necessarily part of a center-surround and therefore not all activities appear correlated. (B) Cumulative PSTH of the whole GrC and GoC populations of the model along a 5 s period. Note that the two PSTH show marked low-frequency oscillations (average 1.8 Hz) around their average level of activity. (C) Autocorrelograms of activity in the GrC and GoC populations and crosscorrelogram of the GrC and GoC populations (in this example the inhibition among GoCs is switched off). Note the high level of correlation in all the three cases on the same main frequency of 1.8 Hz.
In this paper a new scaffold modeling strategy is presented, that is used to simulate fundamental functional properties of the cerebellar microcircuit. The cerebellar scaffold includes the canonical neuron types (GrCs, GoCs, PCs, SCs, BCs, DCNCs), each one with specific geometry of dendritic and axonal fields and with specific convergence/divergence ratios for connectivity. The neurons were purposefully simplified into single point models in order to focus on network connectivity before involving more complex neuronal geometries and properties. The circuit functionality was then tested by applying background activity and burst stimuli and evaluating the network responses. In addition to faithfully reproduce a broad range of experimental observations, the cerebellar scaffold shows the emergence of complex spatiotemporal patterns of activity similar to those observed in vivo and eventually predicts the critical role of local connectome for network functionality.
The Scaffold Design
The scaffold includes two modules: cell placement and connectivity. The first module placed neurons in their corresponding layers according to density values derived from literature. Cell placement exploited a bounded self-avoiding random walk algorithm, except for PCs, which required a placement rule accounting for their regular disposition and quasi-planar non-intersecting dendritic trees. The second module generated microcircuit connectivity by defining the pre- and post-synaptic neurons among those intersecting their dendritic and axonal fields and then establishing the corresponding number of synapses through specific connection probabilities. Geometrical constraints and divergence/convergence ratios derived from literature played a critical role to implement the microcircuit connectome. The distributions of soma distances, both for cell positioning and connectivity, were assessed and provided an internal validation for the network construction processes. The cerebellar scaffold was then implemented using LIF neuron models, whose parameters were tuned to approximate the basal firing and input-output relationships of cerebellar neurons. Finally, functional simulations required the scaffold to be uploaded into a neuro-simulator, either pyNEST or pyNEURON, that worked equivalently well for this purpose.
The scaffold modeling strategy used for the cerebellum microcircuit differs from that used for the cortical microcolumn mostly because here the connectivity rules are based on available statistical and geometrical information rather than on single neuron morphologies and touch-detection (Markram et al., 2015). This allows the scaffold to fully exploit the experimental data available in the cerebellar literature despite incomplete availability of detailed morphological reconstructions of cerebellar neurons. By considering that neuron models based on detailed morphological reconstructions are still unavailable in most neuronal circuits, the strategy adopted here for the cerebellum has a large potential for applicability in a variety of different brain microcircuits. It should be noted that our “intersection-connection” rule is formally similar to the “proximity-connection” rule used for touch-detection in Markram et al. (2015). Eventually, the touch detection strategy could be implemented in the scaffold providing a construction alternative, in which connectivity is directly constrained by neuronal morphology. The advantage would be to specifically connect synapses on specific positions of the dendritic tree, fully exploiting non-linear dendritic computations. Also the data positioning rules could be changed, for example by importing cell positions from the Allen Brain Atlas directly (available at https://portal.bluebrain.epfl.ch/) or using network growing algorithms (Setty et al., 2011; Nguyen et al., 2016).
Simulation and Validation of Cerebellar Network Properties
Following background random inputs and punctuate sensory stimulation, the scaffold model predicted a set of relevant network response properties that matched experimental observations. In the granular layer, the GrC and GoC activity in response to random mossy fiber inputs showed loose synchronicity, as observed in vivo (Pellerin and Lamarre, 1997; Hartmann and Bower, 1998). The GrC and GoC activity in response to bursts in mossy fiber bundles revealed a center-surround organization, as reported in vitro (Mapelli and D'Angelo, 2007; Diwakar et al., 2011; Gandolfi et al., 2014), which was enhanced by aa synapses (Cesana et al., 2013). At the level of molecular layer, the spatial PC discharge patterns depended on the geometry of SC and BC inhibition (Santamaria et al., 2007), and PC burst-pause discharges were generated (Herzfeld et al., 2015). The local PC response was enhanced by granule cell aa, as anticipated by Bower and Woolston (1983), Walter et al. (2009), and Cesana et al. (2013), supporting the vertical organization of GrC-PC transmission (Rokni et al., 2007).
Interestingly, despite the use of simplified LIF neuron models, the observation of these activity patterns suggests that structural constraints play a critical role in determining local neuronal dynamics. In particular, connectivity allows the emergence of center-surrounds in the granular layer and spots of PC activity in the molecular layer. There are several aspects that remain to be assessed and will be easily incorporated into more advanced versions of the cerebellar scaffold.
First of all, assessing the role of non-linear neuronal properties, like intrinsic oscillations, resonance bursting and rebounds, requires to incorporate into the scaffold realistic ionic-channel based neuronal models. Along with this, dendritic computation needs morphologically detailed neuron models that are currently under construction and testing. These include the PCs model (e.g., De Schutter and Bower, 1994; Masoli et al., 2015; Masoli and D'Angelo, 2017), the GrC model (Masoli et al., 2017), the GoC model (Solinas et al., 2007a,b; Kanichay and Silver, 2008), the SC and BC model (currently under construction), the DCN model (Steuber and Jaeger, 2013). Dynamic synapses (Tsodyks and Markram, 1997; Nieus et al., 2006; Migliore et al., 2015) will likewise be incorporated to introduce synaptic strength modulation mechanisms.
Secondly, the scaffold could be used to evaluate the trade-off between computational efficiency and precision. Therefore, the present LIF single point neurons could be substituted by others (extended generalized LIF, E-GLIF) embedding non-linear firing properties (e.g., Brette and Gerstner, 2005; Geminiani et al., 2018) and accounting for synaptic dendritic location by modifying the transmission weight depending on the distance of synapses from the soma (Rössert et al., 2016) or based on experimental data when available.
Thirdly, fully implementing cerebellar connectivity requires the introduction of models of the inferior cerebellar olive (IO) (Libster and Yarom, 2013; De Gruijl et al., 2014). This will complete the DCN-PC-IO cerebellar circuit, allowing the model to simulate oscillations in the olivo-cerebellar circuit, their impact on PC dendritic calcium signaling and computation, and eventually climbing fiber control of plasticity at parallel fiber synapses (Coesmans et al., 2004).
Finally, the addition of novel connections and cells, like the PC to GrC inhibition (Guo et al., 2016) in the anterior cerebellum, the unipolar-brush cell subcircuit in the flocculo-nodular lobe (Mugnaini and Floris, 1994; Subramaniyam et al., 2014), or the DCN to granular layer connections (Gao et al., 2016) will allow to further expand the simulation of cerebellar processing in different cerebellar modules.
The scaffold model was able to reconstruct the complex geometry and neuronal interactions of the cerebellar microcircuit based on intersection-connection rules. Given its architectural design, that puts in series interchangeable programming modules, the scaffold could now be used to plug-in different network configurations into neuronal simulators like e.g., pyNEST and pyNEURON. Both the cell placement algorithm, the neuron model types and the connectivity rules could be substituted to assess different construction strategies and adapted to available data to probe specific functional hypotheses. For example, the connectivity could be recalculated using realistic neuronal morphologies and touch-detection algorithms (proximity-connection rule), as in the cortical microcolumn model (Markram et al., 2015). We envisage that this scaffold modeling strategy, given its versatility, will also be able to host microcircuits different from cerebellum, thus providing a new tool for neurocomputational investigations. It should be noted that the reconstruction procedure is python-based and can be imported in many different simulation frameworks. For example, translation of the scaffold model into PyNN would facilitate neurorobotic and neuromorphic hardware applications (Davison et al., 2008).
SC, EM, CM, and CC have written the code and run the simulations and taken part to the production of text and figures. CC has coordinated the modeling and simulation activity. ED has directed the work and edited the manuscript.
This research was supported by the HBP Brain Simulation Platform funded from the European Union's Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 720270 (Human Brain Project SGA1) and under the Specific Grant Agreement No. 785907 (Human Brain Project SGA2).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fninf.2019.00037/full#supplementary-material
Video S1. The movie displays the sequential placement of Stellate cells in the volume following the bounded self-avoiding random-walk algorithm.
Aizenman, C. D., Huang, E. J., and Linden, D. J. (2003). Morphological correlates of intrinsic electrical excitability in neurons of the deep cerebellar nuclei. J. Neurophysiol. 89, 1738–1747. doi: 10.1152/jn.01043.2002
Billings, G., Piasini, E., Lorincz, A., Nusser, Z., and Silver, R. A. (2014). Network structure within the cerebellar input layer enables lossless sparse encoding. Neuron 83, 960–974. doi: 10.1016/j.neuron.2014.07.020
Boele, H. J., Koekkoek, S. K., De Zeeuw, C. I., and Ruigrok, T. J. (2013). Axonal sprouting and formation of terminals in the adult cerebellum during associative motor learning. J. Neurosci. 33, 17897–17907. doi: 10.1523/JNEUROSCI.0511-13.2013
Bower, J. M., and Woolston, D. C. (1983). Congruence of spatial organization of tactile projections to granule cell and Purkinje cell layers of cerebellar hemispheres of the albino rat: vertical organization of cerebellar cortex. J. Neurophysiol. 49, 745–766. doi: 10.1152/jn.19220.127.116.115
Brette, R., Rudolph, M., Carnevale, T., Hines, M., Beeman, D., Bower, J. M., et al. (2007). Simulation of networks of spiking neurons: a review of tools and strategies. J. Comput. Neurosci. 23, 349–398. doi: 10.1007/s10827-007-0038-6
Canepari, M., Papageorgiou, G., Corrie, J. E., Watkins, C., and Ogden, D. (2001). The conductance underlying the parallel fibre slow EPSP in rat cerebellar Purkinje neurones studied with photolytic release of L-glutamate. J. Physiol. 533, 765–772. doi: 10.1111/j.1469-7793.2001.00765.x
Cao, Y., Maran, S. K., Dhamala, M., Jaeger, D., and Heck, D. H. (2012). Behavior-related pauses in simple-spike activity of mouse Purkinje cells are linked to spike rate modulation. J. Neurosci. 32, 8678–8685. doi: 10.1523/JNEUROSCI.4969-11.2012
Cesana, E., Pietrajtis, K., Bidoret, C., Isope, P., D'Angelo, E., Dieudonne, S., et al. (2013). Granule cell ascending axon excitatory synapses onto Golgi cells implement a potent feedback circuit in the cerebellar granular layer. J. Neurosci. 33, 12430–12446. doi: 10.1523/JNEUROSCI.4897-11.2013
Chu, C. P., Bing, Y. H., Liu, H., and Qiu, D. L. (2012). Roles of molecular layer interneurons in sensory information processing in mouse cerebellar cortex Crus II in vivo. PLoS ONE 7:e37031. doi: 10.1371/journal.pone.0037031
Coesmans, M., Weber, J. T., De Zeeuw, C. I., and Hansel, C. (2004). Bidirectional parallel fiber plasticity in the cerebellum under climbing fiber control. Neuron 44, 691–700. doi: 10.1016/j.neuron.2004.10.031
D'Angelo, E., Antonietti, A., Casali, S., Casellato, C., Garrido, J. A., Luque, N. R., et al. (2016). Modeling the cerebellar microcircuit: new strategies for a long-standing issue. Front. Cell. Neurosci. 10:176. doi: 10.3389/fncel.2016.00176
D'Angelo, E., Nieus, T., Maffei, A., Armano, S., Rossi, P., Taglietti, V., et al. (2001). Theta-frequency bursting and resonance in cerebellar granule cells: experimental evidence and modeling of a slow k+-dependent mechanism. J. Neurosci. 21, 759–770. doi: 10.1523/JNEUROSCI.21-03-00759.2001
D'Angelo, E., Solinas, S., Mapelli, J., Gandolfi, D., Mapelli, L., and Prestori, F. (2013). The cerebellar Golgi cell and spatiotemporal organization of granular layer activity. Front. Neural Circuits 7:93. doi: 10.3389/fncir.2013.00093
D'Angelo, E., and Wheeler-Kingshott, C. (2017). “Modelling the brain: elementary components to explain ensemble functions,” in Rivista Del Nuovo Cimento, ed L. Cifarelli (Bologna: Società Italiana di Fisica, Italian Physical Society), 36.
Davison, A. P., Bruderle, D., Eppler, J., Kremkow, J., Muller, E., Pecevski, D., et al. (2008). PyNN: a common interface for neuronal network simulators. Front. Neuroinform. 2:11. doi: 10.3389/neuro.11.011.2008
De Gruijl, J. R., Sokol, P. A., Negrello, M., and De Zeeuw, C. I. (2014). Modulation of electrotonic coupling in the inferior olive by inhibitory and excitatory inputs: integration in the glomerulus. Neuron 81, 1215–1217. doi: 10.1016/j.neuron.2014.03.009
De Schutter, E., and Bower, J. M. (1994). An active membrane model of the cerebellar Purkinje cell. I. Simulation of current clamps in slice. J. Neurophysiol. 71, 375–400. doi: 10.1152/jn.1918.104.22.1685
Destexhe, A., Bal, T., McCormick, D. A., and Sejnowski, T. J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. 76, 2049–2070. doi: 10.1152/jn.1922.214.171.1249
Dieudonne, S. (1998). Submillisecond kinetics and low efficacy of parallel fibre-Golgi cell synaptic currents in the rat cerebellum. J. Physiol. 510(Pt. 3), 845–866. doi: 10.1111/j.1469-7793.1998.845bj.x
Diwakar, S., Lombardo, P., Solinas, S., Naldi, G., and D'Angelo, E. (2011). Local field potential modeling predicts dense activation in cerebellar granule cells clusters under LTP and LTD control. PLoS ONE 6:e21928. doi: 10.1371/journal.pone.0021928
Gandolfi, D., Pozzi, P., Tognolina, M., Chirico, G., Mapelli, J., and D'Angelo, E. (2014). The spatiotemporal organization of cerebellar network activity resolved by two-photon imaging of multiple single neurons. Front. Cell. Neurosci. 8:92. doi: 10.3389/fncel.2014.00092
Gao, Z., Proietti-Onori, M., Lin, Z., Ten Brinke, M. M., Boele, H. J., Potters, J. W., et al. (2016). Excitatory cerebellar nucleocortical circuit provides internal amplification during associative conditioning. Neuron 89, 645–657. doi: 10.1016/j.neuron.2016.01.008
Geminiani, A., Casellato, C., Locatelli, F., Prestori, F., Pedrocchi, A., and D'Angelo, E. (2018). Complex dynamics in simplified neuronal models: reproducing golgi cell electroresponsiveness. Front. Neuroinform. 12:88. doi: 10.3389/fninf.2018.00088
Gundappa-Sulur, G., De Schutter, E., and Bower, J. M. (1999). Ascending granule cell axon: an important component of cerebellar cortical circuitry. J. Comp. Neurol. 408, 580–596. doi: 10.1002/(SICI)1096-9861(19990614)408:4<580::AID-CNE11>3.0.CO;2-O
Guo, C., Witter, L., Rudolph, S., Elliott, H. L., Ennis, K. A., and Regehr, W. G. (2016). Purkinje cells directly inhibit granule cells in specialized regions of the cerebellar cortex. Neuron 91, 1330–1341. doi: 10.1016/j.neuron.2016.08.011
Heine, S. A., Highstein, S. M., and Blazquez, P. M. (2010). Golgi cells operate as state-specific temporal filters at the input stage of the cerebellar cortex. J. Neurosci. 30, 17004–17014. doi: 10.1523/JNEUROSCI.3513-10.2010
Jaeger, D., and Bower, J. M. (1994). Prolonged responses in rat cerebellar Purkinje cells following activation of the granule cell layer: an intracellular in vitro and in vivo investigation. Exp. Brain Res. 100, 200–214. doi: 10.1007/BF00227191
Jorntell, H., Bengtsson, F., Schonewille, M., and De Zeeuw, C. I. (2010). Cerebellar molecular layer interneurons—computational properties and roles in learning. Trends Neurosci. 33, 524–532. doi: 10.1016/j.tins.2010.08.004
Kanichay, R. T., and Silver, R. A. (2008). Synaptic and cellular properties of the feedforward inhibitory circuit within the input layer of the cerebellar cortex. J. Neurosci. 28, 8955–8967. doi: 10.1523/JNEUROSCI.5469-07.2008
Korbo, L., Andersen, B. B., Ladefoged, O., and Møller, A. (1993). Total numbers of various cell types in rat cerebellar cortex estimated using an unbiased stereological method. Brain Res. 609, 262–268. doi: 10.1016/0006-8993(93)90881-M
Lennon, W., Hecht-Nielsen, R., and Yamazaki, T. (2014). A spiking network model of cerebellar Purkinje cells and molecular layer interneurons exhibiting irregular firing. Front. Comput. Neurosci. 8:157. doi: 10.3389/fncom.2014.00157
Maex, R., and De Schutter, E. (1998). Synchronization of golgi and granule cell firing in a detailed network model of the cerebellar granule cell layer. J. Neurophysiol. 80, 2521–2537. doi: 10.1152/jn.19126.96.36.1991
Maex, R., Vos, B. P., and De Schutter, E. (2000). Weak common parallel fibre synapses explain the loose synchrony observed between rat cerebellar golgi cells. J. Physiol. 523(Pt. 1), 175–192. doi: 10.1111/j.1469-7793.2000.t01-1-00175.x
Markram, H., Muller, E., Ramaswamy, S., Hill, S. L., Segev, I., Schürmann, F., et al. (2015). Reconstruction and simulation of neocortical microcircuitry. Cell 163, 456–492. doi: 10.1016/j.cell.2015.09.029
Masoli, S., and D'Angelo, E. (2017). Synaptic activation of a detailed purkinje cell model predicts voltage-dependent control of burst-pause responses in active dendrites. Front. Cell. Neurosci. 11:278. doi: 10.3389/fncel.2017.00278
Masoli, S., Rizza, M. F., Sgritta, M., Van Geit, W., Schurmann, F., and D'Angelo, E. (2017). Single neuron optimization as a basis for accurate biophysical modeling: the case of cerebellar granule cells. Front. Cell. Neurosci. 11:71. doi: 10.3389/fncel.2017.00071
Masoli, S., Solinas, S., and D'Angelo, E. (2015). Action potential processing in a detailed Purkinje cell model reveals a critical role for axonal compartmentalization. Front. Cell. Neurosci. 9:47. doi: 10.3389/fncel.2015.00047
Migliore, M., Cavarretta, F., Marasco, A., Tulumello, E., Hines, M. L., and Shepherd, G. M. (2015). Synaptic clusters function as odor operators in the olfactory bulb. Proc. Natl. Acad. Sci. U.S.A. 112, 8499–8504. doi: 10.1073/pnas.1502513112
Nguyen, H., Dayan, P., Pujic, Z., Cooper-White, J., and Goodhill, G. J. (2016). A mathematical model explains saturating axon guidance responses to molecular gradients. Elife 5:e12248. doi: 10.7554/eLife.12248
Nieus, T., Sola, E., Mapelli, J., Saftenku, E., Rossi, P., and D'Angelo, E. (2006). LTP regulates burst initiation and frequency at mossy fiber-granule cell synapses of rat cerebellum: experimental observations and theoretical predictions. J. Neurophysiol. 95, 686–699. doi: 10.1152/jn.00696.2005
Ramakrishnan, K. B., Voges, K., De Propris, L., De Zeeuw, C. I., and D'Angelo, E. (2016). Tactile stimulation evokes long-lasting potentiation of purkinje cell discharge in vivo. Front. Cell. Neurosci. 10:36. doi: 10.3389/fncel.2016.00036
Rancz, E. A., Ishikawa, T., Duguid, I., Chadderton, P., Mahon, S., and Hausser, M. (2007). High-fidelity transmission of sensory information by single cerebellar mossy fibre boutons. Nature 450, 1245–1248. doi: 10.1038/nature05995
Roggeri, L., Rivieccio, B., Rossi, P., and D'Angelo, E. (2008). Tactile stimulation evokes long-term synaptic plasticity in the granular layer of cerebellum. J. Neurosci. 28, 6354–6359. doi: 10.1523/JNEUROSCI.5709-07.2008
Santamaria, F., Tripp, P. G., and Bower, J. M. (2007). Feedforward inhibition controls the spread of granule cell-induced Purkinje cell activity in the cerebellar cortex. J. Neurophysiol. 97, 248–263. doi: 10.1152/jn.01098.2005
Setty, Y., Chen, C. C., Secrier, M., Skoblov, N., Kalamatianos, D., and Emmott, S. (2011). How neurons migrate: a dynamic in-silico model of neuronal migration in the developing cortex. BMC Syst. Biol. 5:154. doi: 10.1186/1752-0509-5-154
Sims, R. E., and Hartell, N. A. (2005). Differences in transmission properties and susceptibility to long-term depression reveal functional specialization of ascending axon and parallel fiber synapses to Purkinje cells. J. Neurosci. 25, 3246–3257. doi: 10.1523/JNEUROSCI.0073-05.2005
Solinas, S., Forti, L., Cesana, E., Mapelli, J., De Schutter, E., and D'Angelo, E. (2007a). Fast-reset of pacemaking and theta-frequency resonance patterns in cerebellar golgi cells: simulations of their impact in vivo. Front. Cell. Neurosci. 1:4. doi: 10.3389/neuro.03.004.2007
Solinas, S., Forti, L., Cesana, E., Mapelli, J., De Schutter, E., and D'Angelo, E. (2007b). Computational reconstruction of pacemaking and intrinsic electroresponsiveness in cerebellar Golgi cells. Front. Cell. Neurosci. 1:2. doi: 10.3389/neuro.03.002.2007
Solinas, S., Nieus, T., and D'Angelo, E. (2010). A realistic large-scale model of the cerebellum granular layer predicts circuit spatio-temporal filtering properties. Front. Cell. Neurosci. 4:12. doi: 10.3389/fncel.2010.00012
Subramaniyam, S., Solinas, S., Perin, P., Locatelli, F., Masetto, S., and D'Angelo, E. (2014). Computational modeling predicts the ionic mechanism of late-onset responses in unipolar brush cells. Front. Cell. Neurosci. 8:237. doi: 10.3389/fncel.2014.00237
Sudhakar, S. K., Hong, S., Raikov, I., Publio, R., Lang, C., Close, T., et al. (2017). Spatiotemporal network coding of physiological mossy fiber inputs by the cerebellar granular layer. PLoS Comput. Biol. 13:e1005754. doi: 10.1371/journal.pcbi.1005754
Sultan, F. (2001). Distribution of mossy fibre rosettes in the cerebellum of cat and mice: evidence for a parasagittal organization at the single fibre level. Eur. J. Neurosci. 13, 2123–2130. doi: 10.1046/j.0953-816x.2001.01593.x
Sultan, F., and Heck, D. (2003). Detection of sequences in the cerebellar cortex: numerical estimate of the possible number of tidal-wave inducing sequences represented. J. Physiol. Paris 97, 591–600. doi: 10.1016/j.jphysparis.2004.01.016
Sunkin, S. M., Ng, L., Lau, C., Dolbeare, T., Gilbert, T. L., Thompson, C. L., et al. (2013). Allen Brain Atlas: an integrated spatio-temporal portal for exploring the central nervous system. Nucleic Acids Res. 41, D996–D1008. doi: 10.1093/nar/gks1042
Tripathy, S. J., Savitskaya, J., Burton, S. D., Urban, N. N., and Gerkin, R. C. (2014). NeuroElectro: a window to the world's neuron electrophysiology data. Front. Neuroinform. 8:40. doi: 10.3389/fninf.2014.00040
Tsodyks, M. V., and Markram, H. (1997). The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability. Proc. Natl. Acad. Sci. U.S.A. 94, 719–723. doi: 10.1073/pnas.94.2.719
Walter, J. T., Dizon, M. J., and Khodakhah, K. (2009). The functional equivalence of ascending and parallel fiber inputs in cerebellar computation. J. Neurosci. 29, 8462–8473. doi: 10.1523/JNEUROSCI.5718-08.2009
Keywords: cerebellum, computational spiking models, Python, pyNEST, pyNEURON, connectome
Citation: Casali S, Marenzi E, Medini C, Casellato C and D'Angelo E (2019) Reconstruction and Simulation of a Scaffold Model of the Cerebellar Network. Front. Neuroinform. 13:37. doi: 10.3389/fninf.2019.00037
Received: 21 December 2018; Accepted: 29 April 2019;
Published: 15 May 2019.
Edited by:Gaute T. Einevoll, Norwegian University of Life Sciences, Norway
Reviewed by:Michele Giugliano, University of Antwerp, Belgium
Tadashi Yamazaki, University of Electro-Communications, Japan
Copyright © 2019 Casali, Marenzi, Medini, Casellato and D'Angelo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†These authors have contributed equally to this work as co-last authors